---
  title: "Replication File: On the Mounds of Inequality"
output:
  word_document: default
html_notebook: default
pdf_document: default
---
  
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



Section 1: Set working directory, install and load packages, and program functions necessary for the data analysis.



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 

1.  Set working directory.

```{r}
getwd()
setwd(...)
```

2.  Install and load the relevant packages.

```{r}
# install.packages("EnvStats")
# install.packages("readxl")
# install.packages("plyr")
# install.packages("ineq")
# install.packages("writexl")
# install.packages("xtable")
# install.packages("plotrix")
# install.packages("ggplot2")
# install.packages("plotly")
# install.packages("ggpubr")
# install.packages("rmarkdown")
# install.packages("devtools")
# install.packages("ggspatial")
# install.packages("rnaturalearth")
# install.packages("rnaturalearthdata")
# install.packages("sf")
# install.packages("patchwork")

library("EnvStats")
library("plyr")
library("ineq")
library("writexl")
library("xtable")
library("plotrix")
library("ggplot2")
library("plotly")
library("ggpubr")
library("rmarkdown")
library("devtools")
library("readxl")
library("ggspatial")
library("rnaturalearth")
library("rnaturalearthdata")
library("sf")
library("patchwork")
```

3. Load the dataset and rename average dating variable and multiply barrow_volume with the constant factor.

```{r}
data = read.csv("main_dataset.csv", header = TRUE, 
                    stringsAsFactors = FALSE, encoding = "UTF-8", 
                    na.strings = c("NA",""), strip.white = TRUE)

names(x = data)[names(x = data)=="average_dating_org"] = "average_dating"

data$barrow_volume = 2/(3*pi^0.5)*data$barrow_volume
```

4.  Define functions which will be used in the analysis.

4.1 Generalized Entropy Measures

```{r}
GEM = function(x,c){
  
  n = length(x = x) # number of observations in the sample
  
  
  if(c == 0){I = 1/n*sum(log(x = mean(x = x)/x, base = exp(1)))}else{
    if(c == 1){I = 1/n*sum(x/mean(x)*log(x = x/mean(x), base = exp(1)))}else{
      I = 1/n*1/(c^2-c)*sum((x/mean(x))^c-1)
      
    }
  };
  
  return(I)
  
}
```

4.2 Decomposed Generalized Entropy Measures

```{r}
# => Only use clean data with positive values and without NA's


GEM_Decom = function(x, group_var, obs_var, c){
  
  # group_var: the variable which should be classified into subgroups
  # obs_var: the observation variable e.g., income
  
  data = x # data set as a data frame
  names(data)[names(data) == obs_var] = "obs_var" # rename the column of the observation variable to obs_var
  c = c # sensitivity parameter for GEM selection
  n = as.numeric(x = dim(x = data)[1]) # number of observations in the data set
  a = as.numeric(x = which((colnames(data)==group_var)==1)) # number of the column that contains the group_var
  b = as.numeric(x = which((colnames(data)=="obs_var")==1)) # number of the column that contains the obs_var
  d = data[,a] # vector containing the variable which should be classified into groups
  e = as.numeric(x = data[,b]) # vector containing the values of the observation variable
  names = sort(x = unique(x = d)) # alphabetically ordered group names
  z = count(d) # count the number of observations in each subgroup
  z = z[order(z$x),] # alphabetical order to prevent wrong computations
  
  g = as.numeric(x = length(unique(x = d))) # number of groups
  n_s = as.numeric(x = z[,2]) # number of observations in each subgroup
  x_bar = mean(x = e) # sample mean of the observation variable
  group_means = c(1:g) # vector that will be filled with the group means
  
  for (i in 1:g){
    
    group_means[i] = mean(x = as.matrix(x = subset(x = data, data[a] == names[i], select = c(obs_var)))[,1])
    
  }
  
  I_0wg = as.numeric(x = rep(x = 0, times = g)) # just for later usage 
  I_1wg = as.numeric(x = rep(x = 0, times = g)) # just for later usage
  I_cwg = as.numeric(x = rep(x = 0, times = g)) # just for later usage
  
  # computation of the inequality measure in each group
  
  if(c == 0){
    
    for(i in 1:g){
      
      I_0wg[i] = 1/n_s[i]*sum(x = log(x = group_means[i]/as.matrix(x = subset(x = data, data[a] == names[i], select = c(obs_var)))[,1], base = exp(1)))
      
    }
    
  }else{
    
    if(c == 1){
      
      for(i in 1:g){
        
        I_1wg[i] = 1/n_s[i]*sum(x = as.matrix(x = subset(x = data, data[a] == names[i], select = c(obs_var)))[,1]/group_means[i]*log(x = as.matrix(x = subset(x = data, data[a] == names[i], select = c(obs_var)))[,1]/group_means[i], base = exp(1)))
        
      }
      
    }else{
      
      for(i in 1:g){
        
        I_cwg[i] = 1/n_s[i]*1/(c^2-c)*sum(x = (as.matrix(x = subset(x = data, data[a] == names[i], select = c(obs_var)))[,1]/group_means[i])^c-1)
        
      }
      
    }
  }
  
  # computation of between and within inequality 
  
  if(c == 0){
    
    results = c(I_between = sum(x = n_s/n*log(x = x_bar/group_means, base = exp(1))), I_within = sum(x = n_s/n*I_0wg))
    
  }else{
    
    if(c == 1){
      
      results = c(I_between = sum(x = n_s/n*group_means/x_bar*log(x = group_means/x_bar, base = exp(1))), I_within = sum(x = n_s/n*group_means/x_bar*I_1wg))
      
    }else{
      
      results = c(I_between = 1/(c^2-c)*(sum(x = n_s/n*(group_means/x_bar)^c)-1), I_within = sum(x = n_s/n*(group_means/x_bar)^c*I_cwg))
      
    }  
    
  }
  
  
  return(results)
}
```

4.3 Asymptotic Standard Error GEM

```{r}
ASE_GEM = function(x,c){
  
  x_bar = mean(x = x)
  n = length(x = x)
  
  if(c==0){
    
    Z = (x/x_bar)*log(x = x, base = exp(1))
    Z_bar = mean(x = Z)
    
    ASE = (1/n^2*sum(x = (Z-Z_bar)^2))^0.5
    
  }else{
    
    if(c==1){
      
      Z = (x/x_bar)*(log(x = x/x_bar, base = exp(1))-GEM(x = x, c = c)-1)
      Z_bar = mean(x = Z)
      
      ASE = (1/n^2*sum(x = (Z-Z_bar)^2))^0.5
      
    }else{
      
      Z = 1/(c^2-c)*(x/x_bar)^c-c*(x/x_bar)*(GEM(x = x, c = c)+1/(c^2-c))
      Z_bar = mean(x = Z)
      
      ASE = (1/n^2*sum(x = (Z-Z_bar)^2))^0.5
      
    }
    
  }
  
  return(ASE)
  
}
```

4.4 Normalized Gini formula (values between 0 and 1)

```{r}
Gini_norm = function(x){
  
  s = sort(x = x, decreasing = FALSE)
  n = length(x = s)
  
  G_norm = (2/n*(sum(x = s)^(-1)*(sum(x = c(1:n)*s)))-1/n-1)*n/(n-1);
  return(G_norm)
  
}
```

4.5 Asymptotic Standard Error Gini

```{r}
ASE_Gini = function(x){
  
  s = sort(x = x, decreasing = FALSE)
  x_bar = mean(x = s)
  n = length(x = s)
  index = c(1:n)
  Z = c(1:n)
  G = Gini_norm(x = s)
  
  for (i in 1:n) {
    
    Z[i] = -(G+1)*s[i]+(2*index[i]-1)/n*s[i]-2/n*sum(x = s[1:i])
    
  }
  
  Z_bar = mean(x = Z)
  
  ASE = (1/(n*x_bar)^2*sum(x = (Z-Z_bar)^2))^0.5;
  return(ASE)
  
}
```

4.6 Permutation Test for the Gini Index

```{r}
Perm_Gini = function(X,Y,K,B){
  
  
  # X: one sample
  # Y: another sample
  # K: number of permutations
  # B: number of permutations drawn from K for the analysis
  
  # 1. Compute the mean of each sample and determine size of each sample.
  
  
  x_bar = mean(x = X) # sample mean of the sample X
  y_bar = mean(x = Y) # sample mean of the sample Y
  n = length(x = X) # size of sample X
  m = length(x = Y) # size of sample Y
  
  
  # 2. Divide the observations in each sample by the respective sample mean.
  
  
  X_scaled = X/x_bar
  Y_scaled = Y/y_bar
  
  
  # 3. Combine the scaled observations from each sample into one sample. Determine length of Z.
  
  
  Z = c(X_scaled, Y_scaled)
  p = length(x = Z)
  
  
  # 4. Permute this sample K-times.
  
  
  P = matrix(data = 0, nrow = K, ncol = length(x = Z)) # will contain permutations of Z 
  set.seed(1)
  
  for (i in 1:K) {
    
    P[i,] = sample(x = Z, size = length(Z), replace = FALSE)
    
  }
  
  U = unique(x = P) # Contains unique permutations of P. Necessary, to eliminate identical rows in P.
  
  
  # 5. Draw B of the K permutations from 4. Together with the original sample from 3., this gives B+1 samples.
  
  
  set.seed(1)
  index = sample(x = c(1:K), size = B, replace = FALSE)
  
  U = U[index, ]
  
  
  # 6. Compute the test statistic for the original sample from 3.
  
  
  S_0 = (Gini_norm(x = Z[1:n])-Gini_norm(x = Z[(n+1):p]))/(ASE_Gini(x = Z[1:n])^2+ASE_Gini(x =                 Z[(n+1):p])^2)^(0.5)
  
  
  # 7. Compute the test statistic for each of the B permutation samples from 5.
  
  
  S_obs = c(1:B)
  
  for(i in 1:B){
    
    S_obs[i] = (Gini_norm(x = U[i,1:n])-Gini_norm(x = U[i,(n+1):p]))/(ASE_Gini(x =                                    U[i,1:n])^2+ASE_Gini(x = U[i,(n+1):p])^2)^(0.5)
    
  }
  
  
  # 8. Compute the p-value.
  
  
  H = (sum(x = S_obs<=S_0)+1)/(B+1)
  I = (sum(x = S_obs>=S_0)+1)/(B+1)
  
  p_value = 2*min(H,I); # minimal p-value = (1/B)*2
  
  return(p_value)
  
}
```

4.7 Permutation Test for the GEM Indices

```{r}
Perm_GEM = function(X,Y,c,K,B){
  
  
  # X: one sample
  # Y: another sample
  # c: sensitivity parameter
  # K: number of permutations
  # B: number of permutations drawn from K for the analysis
  
  
  # 1. Compute the mean of each sample and determine size of each sample.
  
  
  x_bar = mean(x = X) # sample mean of the sample X
  y_bar = mean(x = Y) # sample mean of the sample Y
  n = length(x = X) # size of sample X
  m = length(x = Y) # size of sample Y
  
  
  # 2. Divide the observations in each sample by the respective sample mean.
  
  
  X_scaled = X/x_bar
  Y_scaled = Y/y_bar
  
  
  # 3. Combine the scaled observations from each sample into one sample. Determine length of Z.
  
  
  Z = c(X_scaled, Y_scaled)
  p = length(x = Z)
  
  
  # 4. Permute this sample K-times.
  
  
  P = matrix(data = 0, nrow = K, ncol = length(x = Z)) # will contain permutations of Z 
  set.seed(1)
  
  for (i in 1:K) {
    
    P[i,] = sample(x = Z, size = length(Z), replace = FALSE)
    
  }
  
  U = unique(x = P) # Contains unique permutations of P. Necessary, to eliminate identical rows in P.
  
  
  # 5. Draw B of the K permutations from 4. Together with the original sample from 3., this gives B+1 samples.
  
  
  set.seed(1)
  index = sample(x = c(1:K), size = B, replace = FALSE)
  
  U = U[index, ]
  
  
  # 6. Compute the test statistic for the original sample from 3.
  
  
  S_0 = (GEM(x = Z[1:n], c = c)-GEM(x = Z[(n+1):p], c = c))/(ASE_GEM(x = Z[1:n], c = c)^2+ASE_GEM(x =          Z[(n+1):p], c = c)^2)^(0.5)
  
  
  # 7. Compute the test statistic for each of the B permutation samples from 5.
  
  
  S_obs = c(1:B)
  
  for(i in 1:B){
    
    S_obs[i] = (GEM(x = U[i,1:n], c = c)-GEM(x = U[i,(n+1):p], c = c))/(ASE_GEM(x = U[i,1:n], c =                     c)^2+ASE_GEM(x = U[i,(n+1):p], c = c)^2)^(0.5)
    
  }
  
  
  # 8. Compute the p-value.
  
  
  H = (sum(x = S_obs<=S_0)+1)/(B+1)
  I = (sum(x = S_obs>=S_0)+1)/(B+1)
  
  p_value = 2*min(H,I); # minimal p-value = (1/B)*2
  
  return(p_value)
  
}
```

4.8 Permutation Test for the GEM Indices - Decomposed Between

```{r}
Perm_GEM_decom_between = function(X,Y,c,group_var,obs_var,K,B){
  
  
  # X: one DATASET
  # Y: another DATASET
  # c: sensitivity parameter
  # group_var: the variable which should be classified into subgroups
  # obs_var: the observation variable e.g., income
  # K: number of permutations
  # B: number of permutations drawn from K for the analysis
  
  
  # 1. Compute the mean of the observation variable in each dataset and determine the number of              observations in each dataset
  
  
  data_X = X # data set as a data frame
  names(data_X)[names(data_X) == obs_var] = "obs_var" # rename the column of the observation variable                                                         to obs_var in dataset X
  names(data_X)[names(data_X) == group_var] = "group_var" # rename the column of the observation                                                                  variable to group_var in dataset X
  
  a_X = as.numeric(x = which((colnames(data_X)=="group_var")==1)) # number of the column that contains                                                                    the group_var
  b_X = as.numeric(x = which((colnames(data_X)=="obs_var")==1)) # number of the column that contains                                                                    the obs_var in dataset X
  e_X = as.numeric(x = data_X[,b_X]) # vector containing the values of the observation variable
  x_bar = mean(x = e_X) # sample mean of the observation variable
  n_X = as.numeric(x = dim(x = data_X)[1]) # number of observations in the dataset X
  
  data_Y = Y # data set as a data frame
  names(data_Y)[names(data_Y) == obs_var] = "obs_var" # rename the column of the observation variable                                                         to obs_var in dataset Y
  names(data_Y)[names(data_Y) == group_var] = "group_var" # rename the column of the observation                                                                  variable to group_var in dataset Y
  a_Y = as.numeric(x = which((colnames(data_Y)=="group_var")==1)) # number of the column that contains                                                                    the group_var
  b_Y = as.numeric(x = which((colnames(data_Y)=="obs_var")==1)) # number of the column that contains                                                                    the obs_var in dataset Y
  e_Y = as.numeric(x = data_Y[,b_Y]) # vector containing the values of the observation variable
  y_bar = mean(x = e_Y) # sample mean of the observation variable
  n_Y = as.numeric(x = dim(x = data_Y)[1]) # number of observations in the dataset Y
  
  
  # 2. Scale the observation variable in each dataset by the respective mean.
  
  
  e_scaled_X = e_X/x_bar
  e_scaled_Y = e_Y/y_bar
  
  data_X[,b_X] = e_scaled_X
  data_Y[,b_Y] = e_scaled_Y
  
  
  # 3. Combine the two datasets. Determine length of Z.
  
  
  Z = rbind(data_X, data_Y)
  p = n_X+n_Y
  
  
  # 4. Reduce the dataset to the columns "group_var" and "obs_var". Permute this new dataset K-times.
  
  Z = as.matrix(Z[,c(a_X,b_X)])
  k = 2*K
  d = seq(from = 1, to = k-1, by = 2);
  
  
  P = matrix(data = 0, nrow = nrow(x = Z), ncol = k) # will contain permutations of Z 
  set.seed(1)
  
  for (i in d) {
    
    P[,c(i,i+1)] = Z[sample(x = nrow(x = Z), size = nrow(x = Z), replace = FALSE),]
    
  }
  
  U = as.data.frame(x = t(x = unique(x = t(x = P)))); # Contains unique permutations of P. Necessary,mto eliminate identical columns in P.
  
  
  # 5. Draw B of the K permutations from 4. Together with the "original dataset" from 3., this gives B+1 datasets.
  
  
  set.seed(1)
  index = sample(x = seq(from = 1, to = k-1, by = 2), size = B, replace = FALSE)
  t = seq(from = 1, to = 2*B, by = 2)+1
  
  U = U[,sort(x = c(index,(index+1)), decreasing = FALSE)]
  U[,t] = sapply(X = U[,t], FUN = as.numeric)
  colnames(x = U) = rep(x = c("group_var","obs_var"), times = B)
  
  # 6. Compute the test statistic for the original sample from 3. Since asymptotic variances are notmavailable for the components, we use the bootstrapped ones instead.
  
  n_boot = 999 # number of bootstraps
  data_X_Boot_between = c(1:n_boot) # empty vector for between components
  data_Y_Boot_between = c(1:n_boot) # empty vector for between components
  set.seed(1)
  
  for(i in 1:length(x = data_X_Boot_between)){
    
    sample_data = data_X[sample(x = nrow(x = data_X), size = nrow(x = data_X), replace = TRUE),]
    
    data_X_Boot_between[i] = as.numeric(x = GEM_Decom(x = sample_data, group_var = "group_var",
                                                      obs_var = "obs_var", c = c)[1])
    
  }
  
  set.seed(1)
  
  for(i in 1:length(x = data_Y_Boot_between)){
    
    sample_data = data_Y[sample(x = nrow(x = data_Y), size = nrow(x = data_Y), replace = TRUE),]
    
    data_Y_Boot_between[i] = as.numeric(x = GEM_Decom(x = sample_data, group_var = "group_var",
                                                      obs_var = "obs_var", c = c)[1])
    
  }
  
  S_0_between = (as.numeric(x = GEM_Decom(x = data_X, group_var = "group_var", obs_var = "obs_var",
                                          c = c)[1])-
                   as.numeric(x = GEM_Decom(x = data_Y, group_var = "group_var", obs_var = "obs_var",
                                            c = c)[1]))/
    (var(x = data_X_Boot_between)+var(x = data_Y_Boot_between))^(0.5)
  
  
  # 7. Compute the test statistic for each of the B permutation samples from 5.
  
  A = 2*B
  C = seq(from = 1, to = A, by = 2)
  S_obs = c(1:A)
  
  B_Var = 75 # number of bootstrap samples used to compute the bootstrapped variances for each of the B test statistics S_obs. Can be adjusted depending on the power of your machine.
  set.seed(1)
  
  for(i in C){
    
    Between_Boot_X = c(1:B_Var)
    Between_Boot_Y = c(1:B_Var)
    Sample_Data_X = U[1:n_X,c(i,(i+1))]
    Sample_Data_Y = U[(n_X+1):(n_X+n_Y),c(i,(i+1))]
    
    for(j in 1:B_Var){
      
      sample_X = Sample_Data_X[sample(x = nrow(Sample_Data_X), size = nrow(Sample_Data_X),
                                      replace = TRUE),]
      sample_Y = Sample_Data_Y[sample(x = nrow(Sample_Data_Y), size = nrow(Sample_Data_Y),
                                      replace = TRUE),]
      
      Between_Boot_X[j] = as.numeric(x = GEM_Decom(x = sample_X, c = c, group_var = "group_var", 
                                                   obs_var = "obs_var")[1])
      Between_Boot_Y[j] = as.numeric(x = GEM_Decom(x = sample_Y, c = c, group_var = "group_var", 
                                                   obs_var = "obs_var")[1])
      
    }
    
    B_Var_X = var(x = Between_Boot_X)
    B_Var_Y = var(x = Between_Boot_Y)
    
    S_obs[i] = (as.numeric(GEM_Decom(x = Sample_Data_X, c = c, group_var = "group_var",
                                     obs_var = "obs_var")[1])+
                  as.numeric(GEM_Decom(x = Sample_Data_Y, c = c, group_var = "group_var",
                                       obs_var = "obs_var")[1]))/(B_Var_X+B_Var_Y)^(0.5)
    
    
    
  }
  
  S_obs = S_obs[!S_obs %in% 1:A]
  
  
  # 8. Compute the p-value.
  
  
  H = (sum(x = S_obs<=S_0_between)+1)/(B+1)
  I = (sum(x = S_obs>=S_0_between)+1)/(B+1)
  
  p_value = 2*min(H,I); # minimal p-value = (1/B)*2
  
  return(p_value)
  
}
```

4.9 Permutation Test for the GEM Indices - Decomposed Within

```{r}
Perm_GEM_decom_within = function(X,Y,c,group_var,obs_var,K,B){
  
  
  # X: one DATASET
  # Y: another DATASET
  # c: sensitivity parameter
  # group_var: the variable which should be classified into subgroups
  # obs_var: the observation variable e.g., income
  # K: number of permutations
  # B: number of permutations drawn from K for the analysis
  
  
  # 1. Compute the mean of the observation variable in each dataset and determine the number of observations in each dataset
  
  
  data_X = X # data set as a data frame
  names(data_X)[names(data_X) == obs_var] = "obs_var" # rename the column of the observation variable                                                         to obs_var in dataset X
  names(data_X)[names(data_X) == group_var] = "group_var" # rename the column of the observation                                                                  variable to group_var in dataset X
  
  a_X = as.numeric(x = which((colnames(data_X)=="group_var")==1)) # number of the column that contains                                                                    the group_var
  b_X = as.numeric(x = which((colnames(data_X)=="obs_var")==1)) # number of the column that contains                                                                    the obs_var in dataset X
  e_X = as.numeric(x = data_X[,b_X]) # vector containing the values of the observation variable
  x_bar = mean(x = e_X) # sample mean of the observation variable
  n_X = as.numeric(x = dim(x = data_X)[1]) # number of observations in the dataset X
  
  data_Y = Y # data set as a data frame
  names(data_Y)[names(data_Y) == obs_var] = "obs_var" # rename the column of the observation variable                                                         to obs_var in dataset Y
  names(data_Y)[names(data_Y) == group_var] = "group_var" # rename the column of the observation                                                                  variable to group_var in dataset Y
  a_Y = as.numeric(x = which((colnames(data_Y)=="group_var")==1)) # number of the column that contains                                                                    the group_var
  b_Y = as.numeric(x = which((colnames(data_Y)=="obs_var")==1)) # number of the column that contains                                                                    the obs_var in dataset Y
  e_Y = as.numeric(x = data_Y[,b_Y]) # vector containing the values of the observation variable
  y_bar = mean(x = e_Y) # sample mean of the observation variable
  n_Y = as.numeric(x = dim(x = data_Y)[1]) # number of observations in the dataset Y
  
  
  # 2. Scale the observation variable in each dataset by the respective mean.
  
  
  e_scaled_X = e_X/x_bar
  e_scaled_Y = e_Y/y_bar
  
  data_X[,b_X] = e_scaled_X
  data_Y[,b_Y] = e_scaled_Y
  
  
  # 3. Combine the two datasets. Determine length of Z.
  
  
  Z = rbind(data_X, data_Y)
  p = n_X+n_Y
  
  
  # 4. Reduce the dataset to the columns "group_var" and "obs_var". Permute this new dataset K-times.
  
  Z = as.matrix(Z[,c(a_X,b_X)])
  k = 2*K
  d = seq(from = 1, to = k-1, by = 2);
  
  
  P = matrix(data = 0, nrow = nrow(x = Z), ncol = k) # will contain permutations of Z 
  set.seed(1)
  
  for (i in d) {
    
    P[,c(i,i+1)] = Z[sample(x = nrow(x = Z), size = nrow(x = Z), replace = FALSE),]
    
  }
  
  U = as.data.frame(x = t(x = unique(x = t(x = P)))); # Contains unique permutations of P. Necessary,                                                         to eliminate identical columns in P.
  
  
  # 5. Draw B of the K permutations from 4. Together with the "original dataset" from 3., this gives B+1 datasets.
  
  
  set.seed(1)
  index = sample(x = seq(from = 1, to = k-1, by = 2), size = B, replace = FALSE)
  t = seq(from = 1, to = 2*B, by = 2)+1
  
  U = U[,sort(x = c(index,(index+1)), decreasing = FALSE)]
  U[,t] = sapply(X = U[,t], FUN = as.numeric)
  colnames(x = U) = rep(x = c("group_var","obs_var"), times = B)
  
  # 6. Compute the test statistic for the original sample from 3. Since asymptotic variances are not available for the components, we use the bootstrapped ones instead.
  
  n_boot = 999 # number of bootstraps
  data_X_Boot_within = c(1:n_boot) # empty vector for between components
  data_Y_Boot_within = c(1:n_boot) # empty vector for between components
  set.seed(1)
  
  for(i in 1:length(x = data_X_Boot_within)){
    
    sample_data = data_X[sample(x = nrow(x = data_X), size = nrow(x = data_X), replace = TRUE),]
    
    data_X_Boot_within[i] = as.numeric(x = GEM_Decom(x = sample_data, group_var = "group_var",
                                                     obs_var = "obs_var", c = c)[2])
    
  }
  
  set.seed(1)
  
  for(i in 1:length(x = data_Y_Boot_within)){
    
    sample_data = data_Y[sample(x = nrow(x = data_Y), size = nrow(x = data_Y), replace = TRUE),]
    
    data_Y_Boot_within[i] = as.numeric(x = GEM_Decom(x = sample_data, group_var = "group_var",
                                                     obs_var = "obs_var", c = c)[2])
    
  }
  
  S_0_within = (as.numeric(x = GEM_Decom(x = data_X, group_var = "group_var", obs_var = "obs_var",
                                         c = c)[2])-
                  as.numeric(x = GEM_Decom(x = data_Y, group_var = "group_var", obs_var = "obs_var",
                                           c = c)[2]))/
    (var(x = data_X_Boot_within)+var(x = data_Y_Boot_within))^(0.5)
  
  
  # 7. Compute the test statistic for each of the B permutation samples from 5.
  
  A = 2*B
  C = seq(from = 1, to = A, by = 2)
  S_obs = c(1:A)
  
  B_Var = 75 # number of bootstrap samples used to compute the bootstrapped variances for each of the B test statistics S_obs. Can be adjusted depending on the power of your machine.
  set.seed(1)
  
  for(i in C){
    
    Within_Boot_X = c(1:B_Var)
    Within_Boot_Y = c(1:B_Var)
    Sample_Data_X = U[1:n_X,c(i,(i+1))]
    Sample_Data_Y = U[(n_X+1):(n_X+n_Y),c(i,(i+1))]
    
    for(j in 1:B_Var){
      
      sample_X = Sample_Data_X[sample(x = nrow(Sample_Data_X), size = nrow(Sample_Data_X),
                                      replace = TRUE),]
      sample_Y = Sample_Data_Y[sample(x = nrow(Sample_Data_Y), size = nrow(Sample_Data_Y),
                                      replace = TRUE),]
      
      Within_Boot_X[j] = as.numeric(x = GEM_Decom(x = sample_X, c = c, group_var = "group_var", 
                                                  obs_var = "obs_var")[2])
      Within_Boot_Y[j] = as.numeric(x = GEM_Decom(x = sample_Y, c = c, group_var = "group_var", 
                                                  obs_var = "obs_var")[2])
      
    }
    
    B_Var_X = var(x = Within_Boot_X)
    B_Var_Y = var(x = Within_Boot_Y)
    
    S_obs[i] = (as.numeric(GEM_Decom(x = Sample_Data_X, c = c, group_var = "group_var",
                                     obs_var = "obs_var")[2])+
                  as.numeric(GEM_Decom(x = Sample_Data_Y, c = c, group_var = "group_var",
                                       obs_var = "obs_var")[2]))/(B_Var_X+B_Var_Y)^(0.5)
    
    
    
  }
  
  S_obs = S_obs[!S_obs %in% 1:A]
  
  
  # 8. Compute the p-value.
  
  
  H = (sum(x = S_obs<=S_0_within)+1)/(B+1)
  I = (sum(x = S_obs>=S_0_within)+1)/(B+1)
  
  p_value = 2*min(H,I); # minimal p-value = (1/B)*2
  
  return(p_value)
  
}
```



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



Section 2: Prepare the dataset for the analysis.



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



1. Delete burial mounds which have a dating length of more than 600 years (pattern robust also to more than 800, 1,000, 1,200, and 1,400 years).

```{r}
# Quantiles of dating intervals before restriction of the dataset.

quantile(data$length_dating)

data = data[-c(which(x = data$length_dating>600)),] # This excludes 1,316 burial mounds from the analysis.

# Quantiles of dating intervals after restriction of the dataset.

quantile(data$length_dating)
```

2. Determine if there are burial mounds which have no barrow form.

```{r}
sum(is.na(data$barrow_form))
```

3. Delete burial mounds which have no area.

```{r}
data[is.na(data)]= 0

data = data[!(data$barrow_area == 0),]

# No loss of data.
```

4. Export the final dataset to Excel.

```{r}
write_xlsx(x = data,
           path = "~/Library/Mobile Documents/com~apple~CloudDocs/CAU Kiel/PhD/Paper_Prehistory_Inequality/Submission_Nature_Humanities/Replication Files/Excel_Files/data.xlsx")
```



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



Section 3: Time Series - 200 years



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



1. Gini Index - Construct the data

```{r}
data_time_200 = data

data_time_200["time_series_dating"] = 1

data_time_200 = data_time_200[-c(which(data_time_200$average_dating>=3800),which(data_time_200$average_dating<0)),]

# 200 year periods for dataset

for(i in 1:dim(x = data_time_200)[1]){
  
    data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<=3800 & data_time_200[i,"average_dating"]>=3600){"3800_3600"}else{data_time_200[i,"time_series_dating"]}
  
    data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<3600 & data_time_200[i,"average_dating"]>=3400){"3600_3400"}else{data_time_200[i,"time_series_dating"]}
  
    data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<3400 & data_time_200[i,"average_dating"]>=3200){"3400_3200"}else{data_time_200[i,"time_series_dating"]}
  
    data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<3200 & data_time_200[i,"average_dating"]>=3000){"3200_3000"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<3000 & data_time_200[i,"average_dating"]>=2800){"3000_2800"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<2800 & data_time_200[i,"average_dating"]>=2600){"2800_2600"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<2600 & data_time_200[i,"average_dating"]>=2400){"2600_2400"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<2400 & data_time_200[i,"average_dating"]>=2200){"2400_2200"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<2200 & data_time_200[i,"average_dating"]>=2000){"2200_2000"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<2000 & data_time_200[i,"average_dating"]>=1800){"2000_1800"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<1800 & data_time_200[i,"average_dating"]>=1600){"1800_1600"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<1600 & data_time_200[i,"average_dating"]>=1400){"1600_1400"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<1400 & data_time_200[i,"average_dating"]>=1200){"1400_1200"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<1200 & data_time_200[i,"average_dating"]>=1000){"1200_1000"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<1000 & data_time_200[i,"average_dating"]>=800){"1000_800"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<800 & data_time_200[i,"average_dating"]>=600){"800_600"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<600 & data_time_200[i,"average_dating"]>=400){"600_400"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<400 & data_time_200[i,"average_dating"]>=200){"400_200"}else{data_time_200[i,"time_series_dating"]}
  
  data_time_200[i,"time_series_dating"] = if(data_time_200[i,"average_dating"]<200 & data_time_200[i,"average_dating"]>=0){"200_0"}else{data_time_200[i,"time_series_dating"]}
  
}

write_xlsx(x = data_time_200,
           path = "~/Library/Mobile Documents/com~apple~CloudDocs/CAU Kiel/PhD/Paper_Prehistory_Inequality/Submission_Nature_Humanities/Replication Files/Excel_Files/data_200_years.xlsx")
```

2. Subset the data with respect to the time intervals.

```{r}
list2env(setNames(split(data_time_200, data_time_200$time_series_dating), 
                  paste0("data_", names(split(data_time_200,data_time_200$time_series_dating)))), envir = .GlobalEnv)
```

3. Plot all burial mounds; irrespective of the time interval.

```{r}
geo_region = ne_countries(scale = "medium", returnclass = "sf")
class(geo_region)

single_map = ggplot(data = geo_region) +
  geom_sf() +
  geom_point(data = data, mapping = aes(x = X, y = Y), size = 1, color = "forestgreen") +
  coord_sf(xlim = c(0, 25), ylim = c(46, 59), expand = FALSE)

pdf(file = "path.pdf",
    width = 34, height = 16)

single_map

dev.off()
```

4. Plot a map for each time interval, merge them in one graph.

```{r}
list_datasets = list(data_200_0,data_400_200,data_600_400,data_800_600,data_1000_800,data_1200_1000,data_1400_1200,data_1600_1400,
                     data_1800_1600,data_2000_1800,data_2200_2000,data_2400_2200,data_2600_2400,data_2800_2600,data_3000_2800,
                     data_3200_3000,data_3400_3200,data_3600_3400,data_3800_3600)

# Define a function that takes a dataset as an argument and returns a ggplot object

map_plot = function(dataset, title){
  
  ggplot(data = geo_region) +
    geom_sf() +
    geom_point(data = dataset, mapping = aes(x = X, y = Y), size = 2, color = "forestgreen") +
    coord_sf(xlim = c(0,25), ylim = c(46, 59), expand = FALSE) +
    ggtitle(sprintf("%s BCE to %s BCE", title[2], title[1])) +
    theme(plot.title = element_text(hjust = 0.5, size = 32),
          axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 26),
          axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 26),
          axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 26),
          axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 26),
          axis.text = element_text(color = "black", size = 26))
  
}

# Define a list of title strings for each dataset

titles = list(c("200","0"),c("400","200"),c("600","400"),c("800","600"),c("1000","800"),c("1200","1000"),c("1400","1200"),
              c("1600","1400"),c("1800","1600"),c("2000","1800"),c("2200","2000"),c("2400","2200"),c("2600","2400"),
              c("2800","2600"),c("3000","2800"),c("3200","3000"),c("3400","3200"),c("3600","3400"),c("3800","3600"))

# Use lapply to create the plots

plots = mapply(map_plot, dataset = list_datasets, title = titles, SIMPLIFY = FALSE)

# Combine the plots into a single plot

combined_maps = plots[[1]]

for (i in 2:length(plots)){
  
  combined_maps = combined_maps + plots[[i]]
  
}

# Display the combined plot

pdf(file = "path.pdf",
    width = 45, height = 30)

combined_maps

dev.off()
```

5. Gini Index - Set up table necessary for the plot and compute the Gini index for each interval.

```{r}
table(data_time_200$time_series_dating)

absolute_periods_200 = c("3800_3600","3600_3400","3400_3200","3200_3000","3000_2800","2800_2600","2600_2400","2400_2200","2200_2000","2000_1800",
                         "1800_1600","1600_1400","1400_1200","1200_1000","1000_800","800_600", "600_400","400_200","200_0")

Gini_time_series_200 = rep(x = 0, times = length(x = absolute_periods_200))

table_time_series_200 = data.frame(absolute_periods_200, Gini_time_series_200)

table_time_series_200["average_date"] = c(3700,3500,3300,3100,2900,2700,2500,2300,2100,1900,1700,1500,1300,1100,900,700,500,300,100)

table_time_series_200["CI_95_G_l"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200["CI_95_G_u"] = rep(1, times = length(x = absolute_periods_200)) 
table_time_series_200["Gini_BSE"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200$n = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200[i,7] = c(length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])]))
  
}

table_time_series_200["sorting"] = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200[i,2] = Gini_norm(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])])
  
}
```

6. Gini Index - Compute the bootstrapped standard errors of the Gini index.

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

Gini_3800_3600 = rep(x = 0, times = B)
Gini_3600_3400 = rep(x = 0, times = B)
Gini_3400_3200 = rep(x = 0, times = B)
Gini_3200_3000 = rep(x = 0, times = B)
Gini_3000_2800 = rep(x = 0, times = B)
Gini_2800_2600 = rep(x = 0, times = B)
Gini_2600_2400 = rep(x = 0, times = B)
Gini_2400_2200 = rep(x = 0, times = B)
Gini_2200_2000 = rep(x = 0, times = B)
Gini_2000_1800 = rep(x = 0, times = B)
Gini_1800_1600 = rep(x = 0, times = B)
Gini_1600_1400 = rep(x = 0, times = B)
Gini_1400_1200 = rep(x = 0, times = B)
Gini_1200_1000 = rep(x = 0, times = B)
Gini_1000_800 = rep(x = 0, times = B)
Gini_800_600 = rep(x = 0, times = B)
Gini_600_400 = rep(x = 0, times = B)
Gini_400_200 = rep(x = 0, times = B)
Gini_200_0 = rep(x = 0, times = B)

s_1 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[1])]
s_2 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[2])]
s_3 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[3])]
s_4 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[4])]
s_5 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[5])]
s_6 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[6])]
s_7 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[7])]
s_8 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[8])]
s_9 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[9])]
s_10 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[10])]
s_11 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[11])]
s_12 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[12])]
s_13 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[13])]
s_14 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[14])]
s_15 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[15])]
s_16 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[16])]
s_17 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[17])]
s_18 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[18])]
s_19 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[19])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3800_3600 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3600_3400 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3400_3200 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3200_3000 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_3000_2800 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2800_2600 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2600_2400 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2400_2200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_2200_2000 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_2000_1800 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1800_1600 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1600_1400 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_1400_1200 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_1200_1000 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_1000_800 = sample(x = s_15, size = length(x = s_15), replace = T)
  Sample_800_600 = sample(x = s_16, size = length(x = s_16), replace = T)
  Sample_600_400 = sample(x = s_17, size = length(x = s_17), replace = T)
  Sample_400_200 = sample(x = s_18, size = length(x = s_18), replace = T)
  Sample_200_0 = sample(x = s_19, size = length(x = s_19), replace = T)
  
  Gini_3800_3600[i] = Gini_norm(x = Sample_3800_3600)
  Gini_3600_3400[i] = Gini_norm(x = Sample_3600_3400)
  Gini_3400_3200[i] = Gini_norm(x = Sample_3400_3200)
  Gini_3200_3000[i] = Gini_norm(x = Sample_3200_3000)
  Gini_3000_2800[i] = Gini_norm(x = Sample_3000_2800)
  Gini_2800_2600[i] = Gini_norm(x = Sample_2800_2600)
  Gini_2600_2400[i] = Gini_norm(x = Sample_2600_2400)
  Gini_2400_2200[i] = Gini_norm(x = Sample_2400_2200)
  Gini_2200_2000[i] = Gini_norm(x = Sample_2200_2000)
  Gini_2000_1800[i] = Gini_norm(x = Sample_2000_1800)
  Gini_1800_1600[i] = Gini_norm(x = Sample_1800_1600)
  Gini_1600_1400[i] = Gini_norm(x = Sample_1600_1400)
  Gini_1400_1200[i] = Gini_norm(x = Sample_1400_1200)
  Gini_1200_1000[i] = Gini_norm(x = Sample_1200_1000)
  Gini_1000_800[i] = Gini_norm(x = Sample_1000_800)
  Gini_800_600[i] = Gini_norm(x = Sample_800_600)
  Gini_600_400[i] = Gini_norm(x = Sample_600_400)
  Gini_400_200[i] = Gini_norm(x = Sample_400_200)
  Gini_200_0[i] = Gini_norm(x = Sample_200_0)
  
}

Ginis_Boot_time_all_200 = rbind(Gini_3800_3600,Gini_3600_3400,Gini_3400_3200,Gini_3200_3000,Gini_3000_2800, Gini_2800_2600, Gini_2600_2400,
                                Gini_2400_2200,Gini_2200_2000, Gini_2000_1800, Gini_1800_1600, Gini_1600_1400,Gini_1400_1200, Gini_1200_1000,
                                Gini_1000_800, Gini_800_600,Gini_600_400, Gini_400_200, Gini_200_0)

for (i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200[i,6] = sd(x = Ginis_Boot_time_all_200[i,])
  
}
```

7. Gini Index - Compute the confidence intervals.

```{r}
table_time_series_200[19,] = 0 # since there are too few observations
Ginis_Boot_time_all_200[19, ] = 0

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200[i,4] = table_time_series_200[i,2]-abs(x = quantile(x = Ginis_Boot_time_all_200[i,], probs = c(0.025))[1])*table_time_series_200[i,6]
  table_time_series_200[i,5] = table_time_series_200[i,2]+abs(x = quantile(x = Ginis_Boot_time_all_200[i,], probs = c(0.975))[1])*table_time_series_200[i,6]
  
}

# Delete rows with too few observations

table_time_series_200 = table_time_series_200[-c(4,5,19),]

table_time_series_200

# Adjust absolute_periods_200 vector

absolute_periods_200_adjusted = c("3800_3600","3600_3400","3400_3200","2800_2600","2600_2400","2400_2200","2200_2000","2000_1800","1800_1600",
                                  "1600_1400","1400_1200","1200_1000","1000_800","800_600", "600_400","400_200")

xtable(table_time_series_200, type = "latex") # for Latex code

write.table(table_time_series_200, file = "summary_Gini.txt", sep = ";", quote = FALSE, row.names = F)
```

8. P-Values 200-Year-Intervals, Gini Index

```{r}
p_values_Intervals_Gini_200 = c(1:15)

for(i in 1:15){
  
  p_values_Intervals_Gini_200[i] = Perm_Gini(X = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i])],
                                             Y = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i+1])],
                                             K = 20000, B = 999)
  
}



# Check statistical significance of the difference between time slice "2400_2200" and "2000_1800"

Perm_Gini(X = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                          absolute_periods_200_adjusted[6])], Y = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[8])], K = 20000, B = 999)

# Check statistical significance of the difference between time slice "1200_1000" and "800_600"

Perm_Gini(X = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                          absolute_periods_200_adjusted[12])], Y = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[14])], K = 20000, B = 999)

p_values_Intervals_Gini_200
```

9. Gini Index - Plot the results.

```{r}
graph_o_200 = table_time_series_200 %>%
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = Gini_time_series_200), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_200, size = 2,
            mapping = aes(x = -average_date, y = Gini_time_series_200), color = "coral2", linetype = "solid") +
  geom_text(data = table_time_series_200, mapping = aes(x = -average_date, y = Gini_time_series_200, 
                                                        label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("G"^"*")) +
  xlab("Years BCE") +
  ggtitle("Normalized Gini Index") +
  scale_x_continuous(limits = c(-3800,-200),
                     breaks = c(-200,-400,-600,-800,-1000,-1200,-1400,-1600,-1800,-2000,-2200,-2400,-2600,-2800,-3000,-3200,-3400,-3600,-3800),
                     labels = c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  annotate(geom = "text", x = -430, y = 0.935, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 0.885, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -830, y = 0.85, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 0.8, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 0.8, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1430, y = 0.69, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1630, y = 0.65, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1830, y = 0.79, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2030, y = 0.79, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2230, y = 0.60, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.51, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2630, y = 0.588, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3030, y = 0.50, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3430, y = 0.64, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3630, y = 0.70, label = "o", hjust = "left", size = 12, colour = "coral2") +
  coord_cartesian(ylim = c(0.3, 1), xlim = c(-3800,-200), clip = "off")

graph_o_200
```

10. GEM Index (lower part) - Set up table necessary for the plot and compute the GEM index (lower part) for each interval.

```{r}
GEM_0_time_series_200 = rep(x = 0, times = length(x = absolute_periods_200))

table_time_series_200_GEM_0 = data.frame(absolute_periods_200, GEM_0_time_series_200)

table_time_series_200_GEM_0["average_date"] = c(3700,3500,3300,3100,2900,2700,2500,2300,2100,1900,1700,1500,1300,1100,900,700,500,300,100)

table_time_series_200_GEM_0["CI_95_GEM_0_l"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200_GEM_0["CI_95_GEM_0_u"] = rep(1, times = length(x = absolute_periods_200)) 
table_time_series_200_GEM_0["GEM_0_BSE"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200_GEM_0$n = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_0[i,7] = c(length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])]))
  
}

table_time_series_200_GEM_0["sorting"] = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_0[i,2] = GEM(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])], c = 0)
  
}
```

11. GEM Index (lower part) - Compute the bootstrapped standard errors of the GEM index (lower part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_0_3800_3600 = rep(x = 0, times = B)
GEM_0_3600_3400 = rep(x = 0, times = B)
GEM_0_3400_3200 = rep(x = 0, times = B)
GEM_0_3200_3000 = rep(x = 0, times = B)
GEM_0_3000_2800 = rep(x = 0, times = B)
GEM_0_2800_2600 = rep(x = 0, times = B)
GEM_0_2600_2400 = rep(x = 0, times = B)
GEM_0_2400_2200 = rep(x = 0, times = B)
GEM_0_2200_2000 = rep(x = 0, times = B)
GEM_0_2000_1800 = rep(x = 0, times = B)
GEM_0_1800_1600 = rep(x = 0, times = B)
GEM_0_1600_1400 = rep(x = 0, times = B)
GEM_0_1400_1200 = rep(x = 0, times = B)
GEM_0_1200_1000 = rep(x = 0, times = B)
GEM_0_1000_800 = rep(x = 0, times = B)
GEM_0_800_600 = rep(x = 0, times = B)
GEM_0_600_400 = rep(x = 0, times = B)
GEM_0_400_200 = rep(x = 0, times = B)
GEM_0_200_0 = rep(x = 0, times = B)

s_1 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[1])]
s_2 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[2])]
s_3 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[3])]
s_4 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[4])]
s_5 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[5])]
s_6 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[6])]
s_7 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[7])]
s_8 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[8])]
s_9 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[9])]
s_10 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[10])]
s_11 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[11])]
s_12 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[12])]
s_13 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[13])]
s_14 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[14])]
s_15 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[15])]
s_16 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[16])]
s_17 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[17])]
s_18 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[18])]
s_19 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[19])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3800_3600 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3600_3400 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3400_3200 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3200_3000 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_3000_2800 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2800_2600 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2600_2400 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2400_2200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_2200_2000 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_2000_1800 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1800_1600 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1600_1400 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_1400_1200 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_1200_1000 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_1000_800 = sample(x = s_15, size = length(x = s_15), replace = T)
  Sample_800_600 = sample(x = s_16, size = length(x = s_16), replace = T)
  Sample_600_400 = sample(x = s_17, size = length(x = s_17), replace = T)
  Sample_400_200 = sample(x = s_18, size = length(x = s_18), replace = T)
  Sample_200_0 = sample(x = s_19, size = length(x = s_19), replace = T)
  
  GEM_0_3800_3600[i] = GEM(x = Sample_3800_3600, c = 0)
  GEM_0_3600_3400[i] = GEM(x = Sample_3600_3400, c = 0)
  GEM_0_3400_3200[i] = GEM(x = Sample_3400_3200, c = 0)
  GEM_0_3200_3000[i] = GEM(x = Sample_3200_3000, c = 0)
  GEM_0_3000_2800[i] = GEM(x = Sample_3000_2800, c = 0)
  GEM_0_2800_2600[i] = GEM(x = Sample_2800_2600, c = 0)
  GEM_0_2600_2400[i] = GEM(x = Sample_2600_2400, c = 0)
  GEM_0_2400_2200[i] = GEM(x = Sample_2400_2200, c = 0)
  GEM_0_2200_2000[i] = GEM(x = Sample_2200_2000, c = 0)
  GEM_0_2000_1800[i] = GEM(x = Sample_2000_1800, c = 0)
  GEM_0_1800_1600[i] = GEM(x = Sample_1800_1600, c = 0)
  GEM_0_1600_1400[i] = GEM(x = Sample_1600_1400, c = 0)
  GEM_0_1400_1200[i] = GEM(x = Sample_1400_1200, c = 0)
  GEM_0_1200_1000[i] = GEM(x = Sample_1200_1000, c = 0)
  GEM_0_1000_800[i] = GEM(x = Sample_1000_800, c = 0)
  GEM_0_800_600[i] = GEM(x = Sample_800_600, c = 0)
  GEM_0_600_400[i] = GEM(x = Sample_600_400, c = 0)
  GEM_0_400_200[i] = GEM(x = Sample_400_200, c = 0)
  GEM_0_200_0[i] = GEM(x = Sample_200_0, c = 0)
  
}

GEM_0_Boot_time_all_200 = rbind(GEM_0_3800_3600,GEM_0_3600_3400,GEM_0_3400_3200,GEM_0_3200_3000,GEM_0_3000_2800,GEM_0_2800_2600,
                                GEM_0_2600_2400,GEM_0_2400_2200,GEM_0_2200_2000,GEM_0_2000_1800,GEM_0_1800_1600,GEM_0_1600_1400,
                                GEM_0_1400_1200,GEM_0_1200_1000,GEM_0_1000_800, GEM_0_800_600, GEM_0_600_400,GEM_0_400_200,GEM_0_200_0)

for (i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_0[i,6] = sd(x = GEM_0_Boot_time_all_200[i,])
  
}
```

12. GEM Index (lower part) - Compute the confidence intervals.

```{r}
table_time_series_200_GEM_0[19,] = 0 # since there are too few observations
GEM_0_Boot_time_all_200[19, ] = 0

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_0[i,4] = table_time_series_200_GEM_0[i,2]-
    abs(x = quantile(x = GEM_0_Boot_time_all_200[i,], probs = c(0.025))[1])*table_time_series_200_GEM_0[i,6]
  table_time_series_200_GEM_0[i,5] = table_time_series_200_GEM_0[i,2]+
    abs(x = quantile(x = GEM_0_Boot_time_all_200[i,], probs = c(0.975))[1])*table_time_series_200_GEM_0[i,6]
  
}

# Delete rows with too few observations

table_time_series_200_GEM_0 = table_time_series_200_GEM_0[-c(4,5,19),]

table_time_series_200_GEM_0

xtable(table_time_series_200_GEM_0, type = "latex") # for Latex code

write.table(table_time_series_200_GEM_0, file = "summary_I_0.txt", sep = ";", quote = FALSE, row.names = F)
```

13. P-Values 200-Year-Intervals, GEM Index (lower part)

```{r}
p_values_Intervals_GEM_0_200 = c(1:15)

for(i in 1:15){
  
  p_values_Intervals_GEM_0_200[i] = Perm_GEM(X = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i])],
                                             Y = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i+1])],
                                             K = 20000, B = 999, c = 0)
  
}

p_values_Intervals_GEM_0_200
```

14. GEM Index (lower part) - Plot the results.

```{r}
graph_p_200 = table_time_series_200_GEM_0 %>%
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_0_time_series_200), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_200_GEM_0, size = 2,
            mapping = aes(x = -average_date, y = GEM_0_time_series_200), color = "coral2", linetype = "solid") +
  geom_text(data = table_time_series_200_GEM_0, mapping = aes(x = -average_date, y = GEM_0_time_series_200, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[0])) +
  xlab("Years BCE") +
  ggtitle("Generalized Entropy Measure (lower part)") +
  scale_x_continuous(limits = c(-3800,-200),
                     breaks = c(-200,-400,-600,-800,-1000,-1200,-1400,-1600,-1800,-2000,-2200,-2400,-2600,-2800,-3000,-3200,-3400,-3600,-3800),
                     labels = c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  annotate(geom = "text", x = -430, y = 2.95, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 2.1, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -830, y = 1.9, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 1.5, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 1.55, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1430, y = 1.1, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1630, y = 0.95, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1830, y = 1.57, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2030, y = 1.53, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2230, y = 0.77, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.64, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2630, y = 0.76, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3030, y = 0.68, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3430, y = 1.3, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3630, y = 1.43, label = "o", hjust = "left", size = 12, colour = "coral2") +
  coord_cartesian(ylim = c(0, 4), xlim = c(-3800,-200), clip = "off")

graph_p_200
```

15. GEM Index (central part) - Set up table necessary for the plot and compute the GEM index (central part) for each interval.

```{r}
GEM_1_time_series_200 = rep(x = 0, times = length(x = absolute_periods_200))

table_time_series_200_GEM_1 = data.frame(absolute_periods_200, GEM_1_time_series_200)

table_time_series_200_GEM_1["average_date"] = c(3700,3500,3300,3100,2900,2700,2500,2300,2100,1900,1700,1500,1300,1100,900,700,500,300,100)

table_time_series_200_GEM_1["CI_95_GEM_1_l"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200_GEM_1["CI_95_GEM_1_u"] = rep(1, times = length(x = absolute_periods_200)) 
table_time_series_200_GEM_1["GEM_1_BSE"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200_GEM_1$n = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_1[i,7] = c(length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])]))
  
}

table_time_series_200_GEM_1["sorting"] = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_1[i,2] = GEM(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])], c = 1)
  
}
```

16. GEM Index (central part) - Compute the bootstrapped standard errors of the GEM index (central part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_1_3800_3600 = rep(x = 0, times = B)
GEM_1_3600_3400 = rep(x = 0, times = B)
GEM_1_3400_3200 = rep(x = 0, times = B)
GEM_1_3200_3000 = rep(x = 0, times = B)
GEM_1_3000_2800 = rep(x = 0, times = B)
GEM_1_2800_2600 = rep(x = 0, times = B)
GEM_1_2600_2400 = rep(x = 0, times = B)
GEM_1_2400_2200 = rep(x = 0, times = B)
GEM_1_2200_2000 = rep(x = 0, times = B)
GEM_1_2000_1800 = rep(x = 0, times = B)
GEM_1_1800_1600 = rep(x = 0, times = B)
GEM_1_1600_1400 = rep(x = 0, times = B)
GEM_1_1400_1200 = rep(x = 0, times = B)
GEM_1_1200_1000 = rep(x = 0, times = B)
GEM_1_1000_800 = rep(x = 0, times = B)
GEM_1_800_600 = rep(x = 0, times = B)
GEM_1_600_400 = rep(x = 0, times = B)
GEM_1_400_200 = rep(x = 0, times = B)
GEM_1_200_0 = rep(x = 0, times = B)

s_1 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[1])]
s_2 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[2])]
s_3 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[3])]
s_4 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[4])]
s_5 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[5])]
s_6 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[6])]
s_7 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[7])]
s_8 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[8])]
s_9 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[9])]
s_10 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[10])]
s_11 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[11])]
s_12 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[12])]
s_13 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[13])]
s_14 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[14])]
s_15 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[15])]
s_16 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[16])]
s_17 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[17])]
s_18 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[18])]
s_19 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[19])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3800_3600 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3600_3400 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3400_3200 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3200_3000 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_3000_2800 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2800_2600 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2600_2400 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2400_2200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_2200_2000 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_2000_1800 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1800_1600 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1600_1400 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_1400_1200 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_1200_1000 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_1000_800 = sample(x = s_15, size = length(x = s_15), replace = T)
  Sample_800_600 = sample(x = s_16, size = length(x = s_16), replace = T)
  Sample_600_400 = sample(x = s_17, size = length(x = s_17), replace = T)
  Sample_400_200 = sample(x = s_18, size = length(x = s_18), replace = T)
  Sample_200_0 = sample(x = s_19, size = length(x = s_19), replace = T)
  
  GEM_1_3800_3600[i] = GEM(x = Sample_3800_3600, c = 1)
  GEM_1_3600_3400[i] = GEM(x = Sample_3600_3400, c = 1)
  GEM_1_3400_3200[i] = GEM(x = Sample_3400_3200, c = 1)
  GEM_1_3200_3000[i] = GEM(x = Sample_3200_3000, c = 1)
  GEM_1_3000_2800[i] = GEM(x = Sample_3000_2800, c = 1)
  GEM_1_2800_2600[i] = GEM(x = Sample_2800_2600, c = 1)
  GEM_1_2600_2400[i] = GEM(x = Sample_2600_2400, c = 1)
  GEM_1_2400_2200[i] = GEM(x = Sample_2400_2200, c = 1)
  GEM_1_2200_2000[i] = GEM(x = Sample_2200_2000, c = 1)
  GEM_1_2000_1800[i] = GEM(x = Sample_2000_1800, c = 1)
  GEM_1_1800_1600[i] = GEM(x = Sample_1800_1600, c = 1)
  GEM_1_1600_1400[i] = GEM(x = Sample_1600_1400, c = 1)
  GEM_1_1400_1200[i] = GEM(x = Sample_1400_1200, c = 1)
  GEM_1_1200_1000[i] = GEM(x = Sample_1200_1000, c = 1)
  GEM_1_1000_800[i] = GEM(x = Sample_1000_800, c = 1)
  GEM_1_800_600[i] = GEM(x = Sample_800_600, c = 1)
  GEM_1_600_400[i] = GEM(x = Sample_600_400, c = 1)
  GEM_1_400_200[i] = GEM(x = Sample_400_200, c = 1)
  GEM_1_200_0[i] = GEM(x = Sample_200_0, c = 1)
  
}

GEM_1_Boot_time_all_200 = rbind(GEM_1_3800_3600,GEM_1_3600_3400,GEM_1_3400_3200,GEM_1_3200_3000,GEM_1_3000_2800,GEM_1_2800_2600,
                                GEM_1_2600_2400,GEM_1_2400_2200,GEM_1_2200_2000,GEM_1_2000_1800,GEM_1_1800_1600,GEM_1_1600_1400,
                                GEM_1_1400_1200,GEM_1_1200_1000,GEM_1_1000_800,GEM_1_800_600,GEM_1_600_400,GEM_1_400_200,GEM_1_200_0)

for (i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_1[i,6] = sd(x = GEM_1_Boot_time_all_200[i,])
  
}
```

17. GEM Index (central part) - Compute the confidence intervals.

```{r}
table_time_series_200_GEM_1[19,] = 0 # since there are too few observations
GEM_1_Boot_time_all_200[19, ] = 0

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_1[i,4] = table_time_series_200_GEM_1[i,2]-
    abs(x = quantile(x = GEM_1_Boot_time_all_200[i,], probs = c(0.025))[1])*table_time_series_200_GEM_1[i,6]
  table_time_series_200_GEM_1[i,5] = table_time_series_200_GEM_1[i,2]+
    abs(x = quantile(x = GEM_1_Boot_time_all_200[i,], probs = c(0.975))[1])*table_time_series_200_GEM_1[i,6]
  
}

# Delete rows with too few observations

table_time_series_200_GEM_1 = table_time_series_200_GEM_1[-c(4,5,19),]

table_time_series_200_GEM_1

xtable(table_time_series_200_GEM_1, type = "latex") # for Latex code

write.table(table_time_series_200_GEM_1, file = "summary_I_1.txt", sep = ";", quote = FALSE, row.names = F)
```

18. P-Values 200-Year-Intervals, GEM (central part)

```{r}
p_values_Intervals_GEM_1_200 = c(1:15)

for(i in 1:15){
  
  p_values_Intervals_GEM_1_200[i] = Perm_GEM(X = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i])],
                                             Y = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i+1])],
                                             K = 20000, B = 999, c = 1)
  
}

p_values_Intervals_GEM_1_200
```

19. GEM Index (central part) - Plot the results.

```{r}
graph_q_200 = table_time_series_200_GEM_1 %>%
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_1_time_series_200), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_200_GEM_1, size = 2,
            mapping = aes(x = -average_date, y = GEM_1_time_series_200), color = "coral2", linetype = "solid") +
  geom_text(data = table_time_series_200_GEM_1, mapping = aes(x = -average_date, y = GEM_1_time_series_200, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[1])) +
  xlab("Years BCE") +
  ggtitle("Generalized Entropy Measure (central part)") +
  scale_x_continuous(limits = c(-3800,-200),
                     breaks = c(-200,-400,-600,-800,-1000,-1200,-1400,-1600,-1800,-2000,-2200,-2400,-2600,-2800,-3000,-3200,-3400,-3600,-3800),
                     labels = c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  annotate(geom = "text", x = -430, y = 3.7, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 2.5, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -830, y = 2.15, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 1.7, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 1.73, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1430, y = 1.02, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1630, y = 0.9, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1830, y = 1.36, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2030, y = 1.32, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2230, y = 0.73, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.6, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2630, y = 0.9, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3030, y = 0.72, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3430, y = 0.8, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3630, y = 0.95, label = "o", hjust = "left", size = 12, colour = "coral2") +
  coord_cartesian(ylim = c(0, 4.4), xlim = c(-3800,-200), clip = "off")

graph_q_200
```

20. GEM Index (upper part) - Set up table necessary for the plot and compute the GEM index (upper part) for each interval.

```{r}
GEM_2_time_series_200 = rep(x = 0, times = length(x = absolute_periods_200))

table_time_series_200_GEM_2 = data.frame(absolute_periods_200, GEM_2_time_series_200)

table_time_series_200_GEM_2["average_date"] = c(3700,3500,3300,3100,2900,2700,2500,2300,2100,1900,1700,1500,1300,1100,900,700,500,300,100)

table_time_series_200_GEM_2["CI_95_GEM_2_l"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200_GEM_2["CI_95_GEM_2_u"] = rep(1, times = length(x = absolute_periods_200)) 
table_time_series_200_GEM_2["GEM_2_BSE"] = rep(1, times = length(x = absolute_periods_200))
table_time_series_200_GEM_2$n = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_2[i,7] = c(length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])]))
  
}

table_time_series_200_GEM_2["sorting"] = c(1:length(x = absolute_periods_200))

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_2[i,2] = GEM(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[i])], c = 2)
  
}
```

21. GEM Index (upper part) - Compute the bootstrapped standard errors of the GEM index (upper part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_2_3800_3600 = rep(x = 0, times = B)
GEM_2_3600_3400 = rep(x = 0, times = B)
GEM_2_3400_3200 = rep(x = 0, times = B)
GEM_2_3200_3000 = rep(x = 0, times = B)
GEM_2_3000_2800 = rep(x = 0, times = B)
GEM_2_2800_2600 = rep(x = 0, times = B)
GEM_2_2600_2400 = rep(x = 0, times = B)
GEM_2_2400_2200 = rep(x = 0, times = B)
GEM_2_2200_2000 = rep(x = 0, times = B)
GEM_2_2000_1800 = rep(x = 0, times = B)
GEM_2_1800_1600 = rep(x = 0, times = B)
GEM_2_1600_1400 = rep(x = 0, times = B)
GEM_2_1400_1200 = rep(x = 0, times = B)
GEM_2_1200_1000 = rep(x = 0, times = B)
GEM_2_1000_800 = rep(x = 0, times = B)
GEM_2_800_600 = rep(x = 0, times = B)
GEM_2_600_400 = rep(x = 0, times = B)
GEM_2_400_200 = rep(x = 0, times = B)
GEM_2_200_0 = rep(x = 0, times = B)

s_1 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[1])]
s_2 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[2])]
s_3 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[3])]
s_4 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[4])]
s_5 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[5])]
s_6 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[6])]
s_7 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[7])]
s_8 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[8])]
s_9 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[9])]
s_10 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[10])]
s_11 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[11])]
s_12 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[12])]
s_13 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[13])]
s_14 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[14])]
s_15 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[15])]
s_16 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[16])]
s_17 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[17])]
s_18 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[18])]
s_19 = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200[19])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3800_3600 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3600_3400 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3400_3200 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3200_3000 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_3000_2800 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2800_2600 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2600_2400 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2400_2200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_2200_2000 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_2000_1800 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1800_1600 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1600_1400 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_1400_1200 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_1200_1000 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_1000_800 = sample(x = s_15, size = length(x = s_15), replace = T)
  Sample_800_600 = sample(x = s_16, size = length(x = s_16), replace = T)
  Sample_600_400 = sample(x = s_17, size = length(x = s_17), replace = T)
  Sample_400_200 = sample(x = s_18, size = length(x = s_18), replace = T)
  Sample_200_0 = sample(x = s_19, size = length(x = s_19), replace = T)
  
  GEM_2_3800_3600[i] = GEM(x = Sample_3800_3600, c = 2)
  GEM_2_3600_3400[i] = GEM(x = Sample_3600_3400, c = 2)
  GEM_2_3400_3200[i] = GEM(x = Sample_3400_3200, c = 2)
  GEM_2_3200_3000[i] = GEM(x = Sample_3200_3000, c = 2)
  GEM_2_3000_2800[i] = GEM(x = Sample_3000_2800, c = 2)
  GEM_2_2800_2600[i] = GEM(x = Sample_2800_2600, c = 2)
  GEM_2_2600_2400[i] = GEM(x = Sample_2600_2400, c = 2)
  GEM_2_2400_2200[i] = GEM(x = Sample_2400_2200, c = 2)
  GEM_2_2200_2000[i] = GEM(x = Sample_2200_2000, c = 2)
  GEM_2_2000_1800[i] = GEM(x = Sample_2000_1800, c = 2)
  GEM_2_1800_1600[i] = GEM(x = Sample_1800_1600, c = 2)
  GEM_2_1600_1400[i] = GEM(x = Sample_1600_1400, c = 2)
  GEM_2_1400_1200[i] = GEM(x = Sample_1400_1200, c = 2)
  GEM_2_1200_1000[i] = GEM(x = Sample_1200_1000, c = 2)
  GEM_2_1000_800[i] = GEM(x = Sample_1000_800, c = 2)
  GEM_2_800_600[i] = GEM(x = Sample_800_600, c = 2)
  GEM_2_600_400[i] = GEM(x = Sample_600_400, c = 2)
  GEM_2_400_200[i] = GEM(x = Sample_400_200, c = 2)
  GEM_2_200_0[i] = GEM(x = Sample_200_0, c = 2)
  
}

GEM_2_Boot_time_all_200 = rbind(GEM_2_3800_3600,GEM_2_3600_3400,GEM_2_3400_3200,GEM_2_3200_3000,GEM_2_3000_2800,GEM_2_2800_2600,
                                GEM_2_2600_2400,GEM_2_2400_2200,GEM_2_2200_2000,GEM_2_2000_1800,GEM_2_1800_1600,GEM_2_1600_1400,
                                GEM_2_1400_1200,GEM_2_1200_1000,GEM_2_1000_800,GEM_2_800_600,GEM_2_600_400,GEM_2_400_200,GEM_2_200_0)

for (i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_2[i,6] = sd(x = GEM_2_Boot_time_all_200[i,])
  
}
```

22. GEM Index (central part) - Compute the confidence intervals.

```{r}
table_time_series_200_GEM_2[19,] = 0 # since there are too few observations
GEM_2_Boot_time_all_200[19, ] = 0

for(i in 1:length(x = absolute_periods_200)){
  
  table_time_series_200_GEM_2[i,4] = table_time_series_200_GEM_2[i,2]-
    abs(x = quantile(x = GEM_1_Boot_time_all_200[i,], probs = c(0.025))[1])*table_time_series_200_GEM_2[i,6]
  table_time_series_200_GEM_2[i,5] = table_time_series_200_GEM_2[i,2]+
    abs(x = quantile(x = GEM_2_Boot_time_all_200[i,], probs = c(0.975))[1])*table_time_series_200_GEM_2[i,6]
  
}

# Delete rows with too few observations

table_time_series_200_GEM_2 = table_time_series_200_GEM_2[-c(4,5,19),]

table_time_series_200_GEM_2

xtable(table_time_series_200_GEM_2, type = "latex") # for Latex code

write.table(table_time_series_200_GEM_2, file = "summary_I_2.txt", sep = ";", quote = FALSE, row.names = F)
```

23. P-Values 200-Year-Intervals, GEM Index (upper part)

```{r}
p_values_Intervals_GEM_2_200 = c(1:15)

for(i in 1:15){
  
  p_values_Intervals_GEM_2_200[i] = Perm_GEM(X = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i])],
                                             Y = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in%                                                 absolute_periods_200_adjusted[i+1])],
                                             K = 20000, B = 999, c = 2)
  
}

p_values_Intervals_GEM_2_200
```

24. GEM Index (upper part) - Plot the results

```{r}
graph_r_200 = table_time_series_200_GEM_2 %>%
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_2_time_series_200), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_200_GEM_2, size = 2,
            mapping = aes(x = -average_date, y = GEM_2_time_series_200), color = "coral2", linetype = "solid") +
  geom_text(data = table_time_series_200_GEM_2, mapping = aes(x = -average_date, y = GEM_2_time_series_200, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[2])) +
  xlab("Years BCE") +
  ggtitle("Generalized Entropy Measure (upper part)") +
  scale_x_continuous(limits = c(-3800,-200),
                     breaks = c(-200,-400,-600,-800,-1000,-1200,-1400,-1600,-1800,-2000,-2200,-2400,-2600,-2800,-3000,-3200,-3400,-3600,-3800),
                     labels = c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  annotate(geom = "text", x = -430, y = 41, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 22, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -830, y = 20, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 11, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1220, y = 11, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1430, y = 4, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1630, y = 3.6, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1830, y = 4.8, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2030, y = 4.7, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2230, y = 3.1, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 3, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2630, y = 4.3, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3030, y = 4, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3430, y = 2.9, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3630, y = 3.2, label = "o", hjust = "left", size = 12, colour = "coral2") +
  coord_cartesian(ylim = c(0, 60), xlim = c(-3800,-200), clip = "off")

graph_r_200
```

25. Plot all graphs in one figure.

```{r}
pdf(file = "path.pdf",
    width = 34, height = 16)

ggarrange(graph_o_200, graph_p_200, graph_q_200, graph_r_200, ncol = 2, nrow = 2)

dev.off()
```

24. Produce descriptive statistics for every 200-year interval.

```{r}
# Set up data frame for summary statistics

stats_3800_3600 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[1])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[1])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[1])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[1])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[1])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[1])]))

stats_3600_3400 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[2])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[2])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[2])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[2])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[2])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[2])]))

stats_3400_3200 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[3])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[3])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[3])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[3])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[3])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[3])]))

stats_2800_2600 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[4])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[4])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[4])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[4])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[4])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[4])]))

stats_2600_2400 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[5])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[5])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[5])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[5])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[5])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[5])]))

stats_2400_2200 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[6])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[6])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[6])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[6])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[6])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[6])]))

stats_2200_2000 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[7])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[7])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[7])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[7])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[7])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[7])]))

stats_2000_1800 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[8])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[8])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[8])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[8])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[8])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[8])]))

stats_1800_1600 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[9])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[9])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[9])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[9])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[9])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[9])]))

stats_1600_1400 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[10])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[10])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[10])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[10])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[10])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[10])]))

stats_1400_1200 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[11])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[11])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[11])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[11])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[11])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[11])]))

stats_1200_1000 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[12])]),
                    max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[12])]),
                    quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[12])],
                             probs = c(0.25,0.5,0.75,0.9)),
                    mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[12])]),
                    sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[12])]),
                    length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[12])]))

stats_1000_800 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[13])]),
                   max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[13])]),
                   quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[13])],
                             probs = c(0.25,0.5,0.75,0.9)),
                   mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[13])]),
                   sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[13])]),
                   length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[13])]))

stats_800_600 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[14])]),
                  max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[14])]),
                  quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[14])],
                             probs = c(0.25,0.5,0.75,0.9)),
                  mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[14])]),
                  sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[14])]),
                  length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[14])]))

stats_600_400 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[15])]),
                  max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[15])]),
                  quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[15])],
                             probs = c(0.25,0.5,0.75,0.9)),
                  mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[15])]),
                  sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[15])]),
                  length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[15])]))

stats_400_200 = c(min(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[16])]),
                  max(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[16])]),
                  quantile(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[16])],
                             probs = c(0.25,0.5,0.75,0.9)),
                  mean(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[16])]),
                  sd(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[16])]),
                  length(x = data_time_200$barrow_volume[which(data_time_200[,"time_series_dating"] %in% absolute_periods_200_adjusted[16])]))


summary_statistics_200_years = round(x =data.frame(rbind(stats_3800_3600,stats_3600_3400,stats_3400_3200,stats_2800_2600,stats_2600_2400,
                                                         stats_2400_2200,stats_2200_2000,stats_2000_1800,stats_1800_1600,stats_1600_1400,
                                                         stats_1400_1200,stats_1200_1000,stats_1000_800,stats_800_600,stats_600_400,
                                                         stats_400_200)),digit = 2)

summary_statistics_200_years$BCE=c(-3700,-3500,-3300,-2700,-2500,-2300,-2100,-1900,-1700,-1500,-1300,-1100,-900,-700,-500,-300)

colnames(summary_statistics_200_years) = c("min","max","q_0.25","q_0.5","q_0.75","q_0.90","mean","std_dev","num_obs","BCE")

xtable(summary_statistics_200_years, type = "latex") # for Latex code

write.table(summary_statistics_200_years, file = "summary_200.txt", sep = ";", quote = FALSE, row.names = F)
```



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



Section 7: A more detailed time series - 250 years



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



1. Gini Index - Construct the data

```{r}
data_time_250 = data

data_time_250["time_series_dating"] = 1

data_time_250 = data_time_250[-c(which(data_time_250$average_dating>=3750),which(data_time_200$average_dating<0)),]

# 250 year periods for dataset

for(i in 1:dim(x = data_time_250)[1]){
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<3750 & data_time_250[i,"average_dating"]>=3500){"3750_3500"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<3500 & data_time_250[i,"average_dating"]>=3250){"3500_3250"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<3250 & data_time_250[i,"average_dating"]>=3000){"3250_3000"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<3000 & data_time_250[i,"average_dating"]>=2750){"3000_2750"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<2750 & data_time_250[i,"average_dating"]>=2500){"2750_2500"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<2500 & data_time_250[i,"average_dating"]>=2250){"2500_2250"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<2250 & data_time_250[i,"average_dating"]>=2000){"2250_2000"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<2000 & data_time_250[i,"average_dating"]>=1750){"2000_1750"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<1750 & data_time_250[i,"average_dating"]>=1500){"1750_1500"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<1500 & data_time_250[i,"average_dating"]>=1250){"1500_1250"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<1250 & data_time_250[i,"average_dating"]>=1000){"1250_1000"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<1000 & data_time_250[i,"average_dating"]>=750){"1000_750"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<750 & data_time_250[i,"average_dating"]>=500){"750_500"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<500 & data_time_250[i,"average_dating"]>=250){"500_250"}else{data_time_250[i,"time_series_dating"]}
  
  data_time_250[i,"time_series_dating"] = if(data_time_250[i,"average_dating"]<250 & data_time_250[i,"average_dating"]>=0){"250_0"}else{data_time_250[i,"time_series_dating"]}
  
}

write_xlsx(x = data_time_250,
           path = "~/Library/Mobile Documents/com~apple~CloudDocs/CAU Kiel/PhD/Paper_Prehistory_Inequality/Submission_Nature_Humanities/Replication Files/Excel_Files/data_250_years.xlsx")
```

2. Gini Index - Set up table necessary for the plot and compute the Gini index for each interval.

```{r}
table(data_time_250$time_series_dating)

absolute_periods_250 = c("3750_3500","3500_3250","3250_3000","3000_2750","2750_2500","2500_2250","2250_2000","2000_1750","1750_1500","1500_1250",
                         "1250_1000","1000_750","750_500","500_250","250_0")

Gini_time_series_250 = rep(x = 0, times = length(x = absolute_periods_250))

table_time_series_250 = data.frame(absolute_periods_250, Gini_time_series_250)

table_time_series_250["average_date"] = c(3625,3375,3125,2875,2625,2375,2125,1875,1625,1375,1125,875,625,375,125)

table_time_series_250["CI_95_G_l"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250["CI_95_G_u"] = rep(1, times = length(x = absolute_periods_250)) 
table_time_series_250["Gini_BSE"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250$n = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250[i,7] = c(length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])]))
  
}

table_time_series_250["sorting"] = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250[i,2] = Gini_norm(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])])
  
}
```

3. Gini Index - Compute the bootstrapped standard errors of the Gini index.

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

Gini_3750_3500 = rep(x = 0, times = B)
Gini_3500_3250 = rep(x = 0, times = B)
Gini_3250_3000 = rep(x = 0, times = B)
Gini_3000_2750 = rep(x = 0, times = B)
Gini_2750_2500 = rep(x = 0, times = B)
Gini_2500_2250 = rep(x = 0, times = B)
Gini_2250_2000 = rep(x = 0, times = B)
Gini_2000_1750 = rep(x = 0, times = B)
Gini_1750_1500 = rep(x = 0, times = B)
Gini_1500_1250 = rep(x = 0, times = B)
Gini_1250_1000 = rep(x = 0, times = B)
Gini_1000_750 = rep(x = 0, times = B)
Gini_750_500 = rep(x = 0, times = B)
Gini_500_250 = rep(x = 0, times = B)
Gini_250_0 = rep(x = 0, times = B)

s_1 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[1])]
s_2 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[2])]
s_3 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[3])]
s_4 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[4])]
s_5 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[5])]
s_6 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[6])]
s_7 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[7])]
s_8 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[8])]
s_9 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[9])]
s_10 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[10])]
s_11 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[11])]
s_12 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[12])]
s_13 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[13])]
s_14 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[14])]
s_15 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[15])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3750_3500 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3500_3250 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3250_3000 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3000_2750 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2750_2500 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2500_2250 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2250_2000 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2000_1750 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1750_1500 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_1500_1250 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1250_1000 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1000_750 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_750_500 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_500_250 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_250_0 = sample(x = s_15, size = length(x = s_15), replace = T)
  
  Gini_3750_3500[i] = Gini_norm(x = Sample_3750_3500)
  Gini_3500_3250[i] = Gini_norm(x = Sample_3500_3250)
  Gini_3250_3000[i] = Gini_norm(x = Sample_3250_3000)
  Gini_3000_2750[i] = Gini_norm(x = Sample_3000_2750)
  Gini_2750_2500[i] = Gini_norm(x = Sample_2750_2500)
  Gini_2500_2250[i] = Gini_norm(x = Sample_2500_2250)
  Gini_2250_2000[i] = Gini_norm(x = Sample_2250_2000)
  Gini_2000_1750[i] = Gini_norm(x = Sample_2000_1750)
  Gini_1750_1500[i] = Gini_norm(x = Sample_1750_1500)
  Gini_1500_1250[i] = Gini_norm(x = Sample_1500_1250)
  Gini_1250_1000[i] = Gini_norm(x = Sample_1250_1000)
  Gini_1000_750[i] = Gini_norm(x = Sample_1000_750)
  Gini_750_500[i] = Gini_norm(x = Sample_750_500)
  Gini_500_250[i] = Gini_norm(x = Sample_500_250)
  Gini_250_0[i] = Gini_norm(x = Sample_250_0)
  
}

Ginis_Boot_time_all_250 = rbind(Gini_3750_3500,Gini_3500_3250,Gini_3250_3000,Gini_3000_2750,Gini_2750_2500,Gini_2500_2250,Gini_2250_2000,
                                Gini_2000_1750,Gini_1750_1500,Gini_1500_1250,Gini_1250_1000,Gini_1000_750,Gini_750_500,Gini_500_250,Gini_250_0)

for (i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250[i,6] = sd(x = Ginis_Boot_time_all_250[i,])
  
}
```

4. Gini Index - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250[i,4] = table_time_series_250[i,2]-abs(x = quantile(x = Ginis_Boot_time_all_250[i,], probs = c(0.025))[1])*table_time_series_250[i,6]
  table_time_series_250[i,5] = table_time_series_250[i,2]+abs(x = quantile(x = Ginis_Boot_time_all_250[i,], probs = c(0.975))[1])*table_time_series_250[i,6]
  
}

table_time_series_250 = table_time_series_250[-c(3,4),] # Too few observations in 3000-2750 BCE

# Adjust absolute_periods_250 vector

absolute_periods_250_adjusted = c("3750_3500","3500_3250","2750_2500","2500_2250","2250_2000","2000_1750","1750_1500","1500_1250","1250_1000",
                                  "1000_750","750_500","500_250","250_0")

xtable(table_time_series_250, type = "latex") # for Latex code

write.table(table_time_series_250, file = "summary_Gini_250.txt", sep = ";", quote = FALSE, row.names = F)
```

5. P-Values 250-Year-Intervals, Gini Index

```{r}
p_values_Intervals_Gini_250 = c(1:12)

for(i in 1:12){
  
  p_values_Intervals_Gini_250[i] = Perm_Gini(X = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i])],
                                             Y = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i+1])], K = 20000, B = 999)
  
}

p_values_Intervals_Gini_250

# Check statistical significance of the difference between time slice "2750_2500" and "2000_1750"

Perm_Gini(X = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                          absolute_periods_250_adjusted[3])], Y = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[6])], K = 20000, B = 999)

# Check statistical significance of the difference between time slice "2000_1750" and "1500_1250"

Perm_Gini(X = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                          absolute_periods_250_adjusted[6])], Y = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[8])], K = 20000, B = 999)
```

6. Gini Index - Plot the results.

```{r}
graph_o_250 = table_time_series_250 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = Gini_time_series_250), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_250, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = Gini_time_series_250), color = "coral2") +
  geom_text(data = table_time_series_250, mapping = aes(x = -average_date, y = Gini_time_series_250, 
                                                        label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("G"^"*")) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3750,0), breaks = c(0,-250,-500,-750,-1000,-1250,-1500,-1750,-2000,-2250,-2500,-2750,-3000,-3250,-3500,-3750),
                     labels = c(0,250,500,750,1000,1250,1500,1750,2000,2250,2500,2750,3000,3250,3500,3750),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Normalized Gini Index") +
  annotate(geom = "text", x = -3530, y = 0.58, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 0.45, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2530, y = 0.52, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2280, y = 0.59, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2030, y = 0.78, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1780, y = 0.8, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 0.68, label = "o", hjust = "left", size = 12, colour = "coral2") +  
  annotate(geom = "text", x = -1280, y = 0.7, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 0.76, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -780, y = 0.87, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -530, y = 0.955, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -280, y = 0.955, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0.3, 1), xlim = c(-3750,0), clip = "off")

graph_o_250
```

7. GEM Index (lower part) - Set up table necessary for the plot and compute the GEM index (lower part) for each interval.

```{r}
GEM_0_time_series_250 = rep(x = 0, times = length(x = absolute_periods_250))

table_time_series_250_GEM_0 = data.frame(absolute_periods_250, GEM_0_time_series_250)

table_time_series_250_GEM_0["average_date"] = c(3625,3375,3125,2875,2625,2375,2125,1875,1625,1375,1125,875,625,375,125)

table_time_series_250_GEM_0["CI_95_GEM_0_l"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250_GEM_0["CI_95_GEM_0_u"] = rep(1, times = length(x = absolute_periods_250)) 
table_time_series_250_GEM_0["GEM_0_BSE"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250_GEM_0$n = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_0[i,7] = c(length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])]))
  
}

table_time_series_250_GEM_0["sorting"] = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_0[i,2] = GEM(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])], c = 0)
  
}
```

8. GEM Index (lower part) - Compute the bootstrapped standard errors of the GEM index (lower part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_0_3750_3500 = rep(x = 0, times = B)
GEM_0_3500_3250 = rep(x = 0, times = B)
GEM_0_3250_3000 = rep(x = 0, times = B)
GEM_0_3000_2750 = rep(x = 0, times = B)
GEM_0_2750_2500 = rep(x = 0, times = B)
GEM_0_2500_2250 = rep(x = 0, times = B)
GEM_0_2250_2000 = rep(x = 0, times = B)
GEM_0_2000_1750 = rep(x = 0, times = B)
GEM_0_1750_1500 = rep(x = 0, times = B)
GEM_0_1500_1250 = rep(x = 0, times = B)
GEM_0_1250_1000 = rep(x = 0, times = B)
GEM_0_1000_750 = rep(x = 0, times = B)
GEM_0_750_500 = rep(x = 0, times = B)
GEM_0_500_250 = rep(x = 0, times = B)
GEM_0_250_0 = rep(x = 0, times = B)

s_1 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[1])]
s_2 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[2])]
s_3 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[3])]
s_4 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[4])]
s_5 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[5])]
s_6 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[6])]
s_7 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[7])]
s_8 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[8])]
s_9 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[9])]
s_10 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[10])]
s_11 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[11])]
s_12 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[12])]
s_13 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[13])]
s_14 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[14])]
s_15 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[15])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3750_3500 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3500_3250 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3250_3000 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3000_2750 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2750_2500 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2500_2250 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2250_2000 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2000_1750 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1750_1500 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_1500_1250 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1250_1000 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1000_750 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_750_500 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_500_250 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_250_0 = sample(x = s_15, size = length(x = s_15), replace = T)
  
  GEM_0_3750_3500[i] = GEM(x = Sample_3750_3500, c = 0)
  GEM_0_3500_3250[i] = GEM(x = Sample_3500_3250, c = 0)
  GEM_0_3250_3000[i] = GEM(x = Sample_3250_3000, c = 0)
  GEM_0_3000_2750[i] = GEM(x = Sample_3000_2750, c = 0)
  GEM_0_2750_2500[i] = GEM(x = Sample_2750_2500, c = 0)
  GEM_0_2500_2250[i] = GEM(x = Sample_2500_2250, c = 0)
  GEM_0_2250_2000[i] = GEM(x = Sample_2250_2000, c = 0)
  GEM_0_2000_1750[i] = GEM(x = Sample_2000_1750, c = 0)
  GEM_0_1750_1500[i] = GEM(x = Sample_1750_1500, c = 0)
  GEM_0_1500_1250[i] = GEM(x = Sample_1500_1250, c = 0)
  GEM_0_1250_1000[i] = GEM(x = Sample_1250_1000, c = 0)
  GEM_0_1000_750[i] = GEM(x = Sample_1000_750, c = 0)
  GEM_0_750_500[i] = GEM(x = Sample_750_500, c = 0)
  GEM_0_500_250[i] = GEM(x = Sample_500_250, c = 0)
  GEM_0_250_0[i] = GEM(x = Sample_250_0, c = 0)
  
}

GEM_0_Boot_time_all_250 = rbind(GEM_0_3750_3500,GEM_0_3500_3250,GEM_0_3250_3000,GEM_0_3000_2750,GEM_0_2750_2500,GEM_0_2500_2250,GEM_0_2250_2000,
                                GEM_0_2000_1750,GEM_0_1750_1500,GEM_0_1500_1250,GEM_0_1250_1000,GEM_0_1000_750,GEM_0_750_500,GEM_0_500_250,
                                GEM_0_250_0)

for (i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_0[i,6] = sd(x = GEM_0_Boot_time_all_250[i,])
  
}
```

9. GEM Index (lower part) - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_0[i,4] = table_time_series_250_GEM_0[i,2]-
    abs(x = quantile(x = GEM_0_Boot_time_all_250[i,], probs = c(0.025))[1])*table_time_series_250_GEM_0[i,6]
  table_time_series_250_GEM_0[i,5] = table_time_series_250_GEM_0[i,2]+
    abs(x = quantile(x = GEM_0_Boot_time_all_250[i,], probs = c(0.975))[1])*table_time_series_250_GEM_0[i,6]
  
}

table_time_series_250_GEM_0 = table_time_series_250_GEM_0[-c(3,4),] # Too few observations in the respective periods

xtable(table_time_series_250_GEM_0, type = "latex") # for Latex code

write.table(table_time_series_250_GEM_0, file = "summary_250_I_0.txt", sep = ";", quote = FALSE, row.names = F)
```

10. P-Values 250-Year-Intervals, GEM Index (lower part)

```{r}
p_values_Intervals_GEM_0_250 = c(1:12)

for(i in 1:12){
  
  p_values_Intervals_GEM_0_250[i] = Perm_GEM(X = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i])],
                                             Y = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i+1])], K = 20000, B = 999, c = 0)
  
}

p_values_Intervals_GEM_0_250
```

11. GEM Index (lower part) - Plot the results.

```{r}
graph_p_250 = table_time_series_250_GEM_0 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_0_time_series_250), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_250_GEM_0, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = GEM_0_time_series_250), color = "coral2") +
  geom_text(data = table_time_series_250_GEM_0, mapping = aes(x = -average_date, y = GEM_0_time_series_250, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[0])) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3750,0), breaks = c(0,-250,-500,-750,-1000,-1250,-1500,-1750,-2000,-2250,-2500,-2750,-3000,-3250,-3500,-3750),
                     labels = c(0,250,500,750,1000,1250,1500,1750,2000,2250,2500,2750,3000,3250,3500,3750),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Generalized Entropy Measure (lower part)") +
  annotate(geom = "text", x = -3530, y = 1, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 0.45, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2530, y = 0.61, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2280, y = 0.7, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2030, y = 1.48, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1780, y = 1.56, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 1.01, label = "o", hjust = "left", size = 12, colour = "coral2") +  
  annotate(geom = "text", x = -1280, y = 1.19, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 1.45, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -780, y = 2, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -530, y = 2.7, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -280, y = 2.6, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0, 3.2), xlim = c(-3750,0), clip = "off")

graph_p_250
```

12. GEM Index (central part) - Set up table necessary for the plot and compute the GEM index (central part) for each interval.

```{r}
GEM_1_time_series_250 = rep(x = 0, times = length(x = absolute_periods_250))

table_time_series_250_GEM_1 = data.frame(absolute_periods_250, GEM_1_time_series_250)

table_time_series_250_GEM_1["average_date"] = c(3625,3375,3125,2875,2625,2375,2125,1875,1625,1375,1125,875,625,375,125)

table_time_series_250_GEM_1["CI_95_GEM_1_l"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250_GEM_1["CI_95_GEM_1_u"] = rep(1, times = length(x = absolute_periods_250)) 
table_time_series_250_GEM_1["GEM_1_BSE"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250_GEM_1$n = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_1[i,7] = c(length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])]))
  
}

table_time_series_250_GEM_1["sorting"] = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_1[i,2] = GEM(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])], c = 1)
  
}
```

13. GEM Index (central part) - Compute the bootstrapped standard errors of the GEM index (central part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_1_3750_3500 = rep(x = 0, times = B)
GEM_1_3500_3250 = rep(x = 0, times = B)
GEM_1_3250_3000 = rep(x = 0, times = B)
GEM_1_3000_2750 = rep(x = 0, times = B)
GEM_1_2750_2500 = rep(x = 0, times = B)
GEM_1_2500_2250 = rep(x = 0, times = B)
GEM_1_2250_2000 = rep(x = 0, times = B)
GEM_1_2000_1750 = rep(x = 0, times = B)
GEM_1_1750_1500 = rep(x = 0, times = B)
GEM_1_1500_1250 = rep(x = 0, times = B)
GEM_1_1250_1000 = rep(x = 0, times = B)
GEM_1_1000_750 = rep(x = 0, times = B)
GEM_1_750_500 = rep(x = 0, times = B)
GEM_1_500_250 = rep(x = 0, times = B)
GEM_1_250_0 = rep(x = 0, times = B)

s_1 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[1])]
s_2 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[2])]
s_3 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[3])]
s_4 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[4])]
s_5 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[5])]
s_6 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[6])]
s_7 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[7])]
s_8 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[8])]
s_9 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[9])]
s_10 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[10])]
s_11 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[11])]
s_12 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[12])]
s_13 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[13])]
s_14 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[14])]
s_15 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[15])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3750_3500 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3500_3250 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3250_3000 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3000_2750 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2750_2500 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2500_2250 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2250_2000 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2000_1750 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1750_1500 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_1500_1250 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1250_1000 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1000_750 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_750_500 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_500_250 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_250_0 = sample(x = s_15, size = length(x = s_15), replace = T)
  
  GEM_1_3750_3500[i] = GEM(x = Sample_3750_3500, c = 1)
  GEM_1_3500_3250[i] = GEM(x = Sample_3500_3250, c = 1)
  GEM_1_3250_3000[i] = GEM(x = Sample_3250_3000, c = 1)
  GEM_1_3000_2750[i] = GEM(x = Sample_3000_2750, c = 1)
  GEM_1_2750_2500[i] = GEM(x = Sample_2750_2500, c = 1)
  GEM_1_2500_2250[i] = GEM(x = Sample_2500_2250, c = 1)
  GEM_1_2250_2000[i] = GEM(x = Sample_2250_2000, c = 1)
  GEM_1_2000_1750[i] = GEM(x = Sample_2000_1750, c = 1)
  GEM_1_1750_1500[i] = GEM(x = Sample_1750_1500, c = 1)
  GEM_1_1500_1250[i] = GEM(x = Sample_1500_1250, c = 1)
  GEM_1_1250_1000[i] = GEM(x = Sample_1250_1000, c = 1)
  GEM_1_1000_750[i] = GEM(x = Sample_1000_750, c = 1)
  GEM_1_750_500[i] = GEM(x = Sample_750_500, c = 1)
  GEM_1_500_250[i] = GEM(x = Sample_500_250, c = 1)
  GEM_1_250_0[i] = GEM(x = Sample_250_0, c = 1)
  
}

GEM_1_Boot_time_all_250 = rbind(GEM_1_3750_3500,GEM_1_3500_3250,GEM_1_3250_3000,GEM_1_3000_2750,GEM_1_2750_2500,GEM_1_2500_2250,GEM_1_2250_2000,
                                GEM_1_2000_1750,GEM_1_1750_1500,GEM_1_1500_1250,GEM_1_1250_1000,GEM_1_1000_750,GEM_1_750_500,GEM_1_500_250,
                                GEM_1_250_0)

for (i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_1[i,6] = sd(x = GEM_1_Boot_time_all_250[i,])
  
}
```

14. GEM Index (central part) - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_1[i,4] = table_time_series_250_GEM_1[i,2]-
    abs(x = quantile(x = GEM_1_Boot_time_all_250[i,], probs = c(0.025))[1])*table_time_series_250_GEM_1[i,6]
  table_time_series_250_GEM_1[i,5] = table_time_series_250_GEM_1[i,2]+
    abs(x = quantile(x = GEM_1_Boot_time_all_250[i,], probs = c(0.975))[1])*table_time_series_250_GEM_1[i,6]
  
}

table_time_series_250_GEM_1 = table_time_series_250_GEM_1[-c(3,4),] # Too few observations in the respective periods

xtable(table_time_series_250_GEM_1, type = "latex") # for Latex code

write.table(table_time_series_250_GEM_1, file = "summary_250_I_1.txt", sep = ";", quote = FALSE, row.names = F)
```

15. P-Values 250-Year-Intervals, GEM Index (central part)

```{r}
p_values_Intervals_GEM_1_250 = c(1:12)

for(i in 1:12){
  
  p_values_Intervals_GEM_1_250[i] = Perm_GEM(X = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i])],
                                             Y = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i+1])], K = 20000, B = 999, c = 1)
  
}

p_values_Intervals_GEM_1_250
```

16. GEM Index (central part) - Plot the results.

```{r}
graph_q_250 = table_time_series_250_GEM_1 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_1_time_series_250), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_250_GEM_1, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = GEM_1_time_series_250), color = "coral2") +
  geom_text(data = table_time_series_250_GEM_1, mapping = aes(x = -average_date, y = GEM_1_time_series_250, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[1])) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3750,0), breaks = c(0,-250,-500,-750,-1000,-1250,-1500,-1750,-2000,-2250,-2500,-2750,-3000,-3250,-3500,-3750),
                     labels = c(0,250,500,750,1000,1250,1500,1750,2000,2250,2500,2750,3000,3250,3500,3750),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Generalized Entropy Measure (central part)") +
  annotate(geom = "text", x = -3530, y = 0.73, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 0.56, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2530, y = 0.65, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2280, y = 0.75, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2030, y = 1.33, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1780, y = 1.4, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 1.01, label = "o", hjust = "left", size = 12, colour = "coral2") +  
  annotate(geom = "text", x = -1280, y = 1.07, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 1.26, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -780, y = 2.3, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -530, y = 3.9, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -280, y = 3.45, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0, 4.5), xlim = c(-3750,0), clip = "off")

graph_q_250
```

17. GEM Index (upper part) - Set up table necessary for the plot and compute the GEM index (upper part) for each interval.

```{r}
GEM_2_time_series_250 = rep(x = 0, times = length(x = absolute_periods_250))

table_time_series_250_GEM_2 = data.frame(absolute_periods_250, GEM_2_time_series_250)

table_time_series_250_GEM_2["average_date"] = c(3625,3375,3125,2875,2625,2375,2125,1875,1625,1375,1125,875,625,375,125)

table_time_series_250_GEM_2["CI_95_GEM_2_l"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250_GEM_2["CI_95_GEM_2_u"] = rep(1, times = length(x = absolute_periods_250)) 
table_time_series_250_GEM_2["GEM_2_BSE"] = rep(1, times = length(x = absolute_periods_250))
table_time_series_250_GEM_2$n = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_2[i,7] = c(length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])]))
  
}

table_time_series_250_GEM_2["sorting"] = c(1:length(x = absolute_periods_250))

for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_2[i,2] = GEM(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[i])], c = 2)
  
}
```

18. GEM Index (upper part) - Compute the bootstrapped standard errors of the GEM index (upper part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_2_3750_3500 = rep(x = 0, times = B)
GEM_2_3500_3250 = rep(x = 0, times = B)
GEM_2_3250_3000 = rep(x = 0, times = B)
GEM_2_3000_2750 = rep(x = 0, times = B)
GEM_2_2750_2500 = rep(x = 0, times = B)
GEM_2_2500_2250 = rep(x = 0, times = B)
GEM_2_2250_2000 = rep(x = 0, times = B)
GEM_2_2000_1750 = rep(x = 0, times = B)
GEM_2_1750_1500 = rep(x = 0, times = B)
GEM_2_1500_1250 = rep(x = 0, times = B)
GEM_2_1250_1000 = rep(x = 0, times = B)
GEM_2_1000_750 = rep(x = 0, times = B)
GEM_2_750_500 = rep(x = 0, times = B)
GEM_2_500_250 = rep(x = 0, times = B)
GEM_2_250_0 = rep(x = 0, times = B)

s_1 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[1])]
s_2 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[2])]
s_3 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[3])]
s_4 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[4])]
s_5 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[5])]
s_6 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[6])]
s_7 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[7])]
s_8 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[8])]
s_9 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[9])]
s_10 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[10])]
s_11 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[11])]
s_12 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[12])]
s_13 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[13])]
s_14 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[14])]
s_15 = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250[15])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3750_3500 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3500_3250 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3250_3000 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_3000_2750 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2750_2500 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2500_2250 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_2250_2000 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_2000_1750 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1750_1500 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_1500_1250 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_1250_1000 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_1000_750 = sample(x = s_12, size = length(x = s_12), replace = T)
  Sample_750_500 = sample(x = s_13, size = length(x = s_13), replace = T)
  Sample_500_250 = sample(x = s_14, size = length(x = s_14), replace = T)
  Sample_250_0 = sample(x = s_15, size = length(x = s_15), replace = T)
  
  GEM_2_3750_3500[i] = GEM(x = Sample_3750_3500, c = 2)
  GEM_2_3500_3250[i] = GEM(x = Sample_3500_3250, c = 2)
  GEM_2_3250_3000[i] = GEM(x = Sample_3250_3000, c = 2)
  GEM_2_3000_2750[i] = GEM(x = Sample_3000_2750, c = 2)
  GEM_2_2750_2500[i] = GEM(x = Sample_2750_2500, c = 2)
  GEM_2_2500_2250[i] = GEM(x = Sample_2500_2250, c = 2)
  GEM_2_2250_2000[i] = GEM(x = Sample_2250_2000, c = 2)
  GEM_2_2000_1750[i] = GEM(x = Sample_2000_1750, c = 2)
  GEM_2_1750_1500[i] = GEM(x = Sample_1750_1500, c = 2)
  GEM_2_1500_1250[i] = GEM(x = Sample_1500_1250, c = 2)
  GEM_2_1250_1000[i] = GEM(x = Sample_1250_1000, c = 2)
  GEM_2_1000_750[i] = GEM(x = Sample_1000_750, c = 2)
  GEM_2_750_500[i] = GEM(x = Sample_750_500, c = 2)
  GEM_2_500_250[i] = GEM(x = Sample_500_250, c = 2)
  GEM_2_250_0[i] = GEM(x = Sample_250_0, c = 2)
  
}

GEM_2_Boot_time_all_250 = rbind(GEM_2_3750_3500,GEM_2_3500_3250,GEM_2_3250_3000,GEM_2_3000_2750,GEM_2_2750_2500,GEM_2_2500_2250,GEM_2_2250_2000,
                                GEM_2_2000_1750,GEM_2_1750_1500,GEM_2_1500_1250,GEM_2_1250_1000,GEM_2_1000_750,GEM_2_750_500,GEM_2_500_250,
                                GEM_2_250_0)

for (i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_2[i,6] = sd(x = GEM_2_Boot_time_all_250[i,])
  
}
```

19. GEM Index (upper part) - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_250)){
  
  table_time_series_250_GEM_2[i,4] = table_time_series_250_GEM_2[i,2]-
    abs(x = quantile(x = GEM_2_Boot_time_all_250[i,], probs = c(0.025))[1])*table_time_series_250_GEM_2[i,6]
  table_time_series_250_GEM_2[i,5] = table_time_series_250_GEM_2[i,2]+
    abs(x = quantile(x = GEM_2_Boot_time_all_250[i,], probs = c(0.975))[1])*table_time_series_250_GEM_2[i,6]
  
}

table_time_series_250_GEM_2 = table_time_series_250_GEM_2[-c(3,4),] # Too few observations in the respective periods

xtable(table_time_series_250_GEM_2, type = "latex") # for Latex code

write.table(table_time_series_250_GEM_2, file = "summary_250_I_2.txt", sep = ";", quote = FALSE, row.names = F)
```

20. P-Values 250-Year-Intervals, GEM Index (upper part)

```{r}
p_values_Intervals_GEM_2_250 = c(1:12)

for(i in 1:12){
  
  p_values_Intervals_GEM_2_250[i] = Perm_GEM(X = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i])],
                                             Y = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in%                                                 absolute_periods_250_adjusted[i+1])], K = 20000, B = 999, c = 2)
  
}

p_values_Intervals_GEM_2_250
```

21. GEM Index (upper part) - Plot the results.

```{r}
graph_r_250 = table_time_series_250_GEM_2 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_2_time_series_250), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_250_GEM_2, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = GEM_2_time_series_250), color = "coral2") +
  geom_text(data = table_time_series_250_GEM_2, mapping = aes(x = -average_date, y = GEM_2_time_series_250, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[2])) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3750,0), breaks = c(0,-250,-500,-750,-1000,-1250,-1500,-1750,-2000,-2250,-2500,-2750,-3000,-3250,-3500,-3750),
                     labels = c(0,250,500,750,1000,1250,1500,1750,2000,2250,2500,2750,3000,3250,3500,3750),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Generalized Entropy Measure (upper part)") +
  annotate(geom = "text", x = -3530, y = 3.5, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 3.6, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2530, y = 4.2, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2280, y = 4.2, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2030, y = 5.2, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1780, y = 5.6, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 4.8, label = "o", hjust = "left", size = 12, colour = "coral2") +  
  annotate(geom = "text", x = -1280, y = 4.8, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 5.2, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -780, y = 23, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -530, y = 62, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -280, y = 51, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0, 78), xlim = c(-3750,0), clip = "off")

graph_r_250
```

22. Plot all graphs in one figure.

```{r}
pdf(file = "path.pdf",
    width = 34, height = 16)

ggarrange(graph_o_250, graph_p_250, graph_q_250, graph_r_250, ncol = 2, nrow = 2)

dev.off()
```

23. Produce descriptive statistics for every 250-year interval.

```{r}
# Set up data frame for summary statistics

stats_3750_3500 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[1])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[1])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[1])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[1])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[1])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[1])]))

stats_3500_3250 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[2])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[2])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[2])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[2])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[2])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[2])]))

stats_2750_2500 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[3])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[3])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[3])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[3])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[3])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[3])]))

stats_2500_2250 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[4])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[4])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[4])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[4])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[4])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[4])]))

stats_2250_2000 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[5])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[5])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[5])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[5])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[5])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[5])]))

stats_2000_1750 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[6])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[6])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[6])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[6])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[6])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[6])]))

stats_1750_1500 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[7])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[7])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[7])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[7])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[7])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[7])]))

stats_1500_1250 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[8])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[8])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[8])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[8])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[8])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[8])]))

stats_1250_1000 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[9])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[9])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[9])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[9])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[9])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[9])]))

stats_1000_750 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[10])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[10])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[10])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[10])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[10])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[10])]))

stats_750_500 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[11])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[11])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[11])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[11])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[11])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[11])]))

stats_500_250 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[12])]),
                    max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[12])]),
                    quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[12])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[12])]),
                    sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[12])]),
                    length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[12])]))

stats_250_0 = c(min(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[13])]),
                   max(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[13])]),
                   quantile(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[13])],
                            probs = c(0.25,0.5,0.75)),
                   mean(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[13])]),
                   sd(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[13])]),
                   length(x = data_time_250$barrow_volume[which(data_time_250[,"time_series_dating"] %in% absolute_periods_250_adjusted[13])]))


summary_statistics_250_years = round(x =data.frame(rbind(stats_3750_3500,stats_3500_3250,stats_2750_2500,stats_2500_2250,stats_2250_2000,
                                                         stats_2000_1750,stats_1750_1500,stats_1500_1250,stats_1250_1000,stats_1000_750,
                                                         stats_750_500,stats_500_250,stats_250_0)),digit = 2)

colnames(summary_statistics_250_years) = c("min","max","q_0.25","q_0.5","q_0.75","mean","std_dev","num_obs")

xtable(summary_statistics_250_years, type = "latex") # for Latex code

write.table(summary_statistics_250_years, file = "summary_250.txt", sep = ";", quote = FALSE, row.names = F)
```



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



Section 8: A more detailed time series - 300 years



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



1. Gini Index - Construct the data

```{r}
data_time_300 = data

data_time_300["time_series_dating"] = 1

# To obtain a proper time series, exclude all burial mounds not in the interval of 3600-0 BCE 

data_time_300 = data_time_300[-c(which(data_time_300$average_dating>3600),which(data_time_300$average_dating<0)),]

# 250 year periods for dataset

for(i in 1:dim(x = data_time_300)[1]){
  
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<3600 & data_time_300[i,"average_dating"]>=3300){"3600_3300"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<3300 & data_time_300[i,"average_dating"]>=3000){"3300_3000"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<3000 & data_time_300[i,"average_dating"]>=2700){"3000_2700"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<2700 & data_time_300[i,"average_dating"]>=2400){"2700_2400"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<2400 & data_time_300[i,"average_dating"]>=2100){"2400_2100"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<2100 & data_time_300[i,"average_dating"]>=1800){"2100_1800"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<1800 & data_time_300[i,"average_dating"]>=1500){"1800_1500"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<1500 & data_time_300[i,"average_dating"]>=1200){"1500_1200"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<1200 & data_time_300[i,"average_dating"]>=900){"1200_900"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<900 & data_time_300[i,"average_dating"]>=600){"900_600"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<600 & data_time_300[i,"average_dating"]>=300){"600_300"}else{data_time_300[i,"time_series_dating"]}
  data_time_300[i,"time_series_dating"] = if(data_time_300[i,"average_dating"]<300 & data_time_300[i,"average_dating"]>=0){"300_0"}else{data_time_300[i,"time_series_dating"]}
  
}

write_xlsx(x = data_time_300,
           path = "~/Library/Mobile Documents/com~apple~CloudDocs/CAU Kiel/PhD/Paper_Prehistory_Inequality/Submission_Nature_Humanities/Replication Files/Excel_files/data_300_years.xlsx")
```

2. Gini Index - Set up table necessary for the plot and compute the Gini index for each interval.

```{r}
table(data_time_300$time_series_dating)

absolute_periods_300 = c("3600_3300","3300_3000","3000_2700","2700_2400","2400_2100","2100_1800","1800_1500","1500_1200","1200_900","900_600",
                         "600_300","300_0")

Gini_time_series_300 = rep(x = 0, times = length(x = absolute_periods_300))

table_time_series_300 = data.frame(absolute_periods_300, Gini_time_series_300)

table_time_series_300["average_date"] = c(3450,3150,2850,2550,2250,1950,1650,1350,1050,750,450,150)

table_time_series_300["CI_95_G_l"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300["CI_95_G_u"] = rep(1, times = length(x = absolute_periods_300)) 
table_time_series_300["Gini_BSE"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300$n = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300[i,7] = c(length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])]))
  
}

table_time_series_300["sorting"] = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300[i,2] = Gini_norm(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])])
  
}
```

3. Gini Index - Compute the bootstrapped standard errors of the Gini index.

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

Gini_3600_3300 = rep(x = 0, times = B)
Gini_3300_3000 = rep(x = 0, times = B)
Gini_3000_2700 = rep(x = 0, times = B)
Gini_2700_2400 = rep(x = 0, times = B)
Gini_2400_2100 = rep(x = 0, times = B)
Gini_2100_1800 = rep(x = 0, times = B)
Gini_1800_1500 = rep(x = 0, times = B)
Gini_1500_1200 = rep(x = 0, times = B)
Gini_1200_900 = rep(x = 0, times = B)
Gini_900_600 = rep(x = 0, times = B)
Gini_600_300 = rep(x = 0, times = B)
Gini_300_0 = rep(x = 0, times = B)

s_1 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]
s_2 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]
s_3 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]
s_4 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]
s_5 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]
s_6 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]
s_7 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]
s_8 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]
s_9 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]
s_10 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]
s_11 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]
s_12 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3600_3300 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3300_3000 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3000_2700 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_2700_2400 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2400_2100 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2100_1800 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_1800_1500 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_1500_1200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1200_900 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_900_600 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_600_300 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_300_0 = sample(x = s_12, size = length(x = s_12), replace = T)
  
  Gini_3600_3300[i] = Gini_norm(x = Sample_3600_3300)
  Gini_3300_3000[i] = Gini_norm(x = Sample_3300_3000)
  Gini_3000_2700[i] = Gini_norm(x = Sample_3000_2700)
  Gini_2700_2400[i] = Gini_norm(x = Sample_2700_2400)
  Gini_2400_2100[i] = Gini_norm(x = Sample_2400_2100)
  Gini_2100_1800[i] = Gini_norm(x = Sample_2100_1800)
  Gini_1800_1500[i] = Gini_norm(x = Sample_1800_1500)
  Gini_1500_1200[i] = Gini_norm(x = Sample_1500_1200)
  Gini_1200_900[i] = Gini_norm(x = Sample_1200_900)
  Gini_900_600[i] = Gini_norm(x = Sample_900_600)
  Gini_600_300[i] = Gini_norm(x = Sample_600_300)
  Gini_300_0[i] = Gini_norm(x = Sample_300_0)
  
}

Ginis_Boot_time_all_300 = rbind(Gini_3600_3300,Gini_3300_3000,Gini_3000_2700,Gini_2700_2400,Gini_2400_2100,Gini_2100_1800,Gini_1800_1500,
                                Gini_1500_1200,Gini_1200_900,Gini_900_600,Gini_600_300,Gini_300_0)

for (i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300[i,6] = sd(x = Ginis_Boot_time_all_300[i,])
  
}
```

4. Gini Index - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300[i,4] = table_time_series_300[i,2]-abs(x = quantile(x = Ginis_Boot_time_all_300[i,], probs = c(0.025))[1])*table_time_series_300[i,6]
  table_time_series_300[i,5] = table_time_series_300[i,2]+abs(x = quantile(x = Ginis_Boot_time_all_300[i,], probs = c(0.975))[1])*table_time_series_300[i,6]
  
}

table_time_series_300

xtable(table_time_series_300, type = "latex") # for Latex code

write.table(table_time_series_300, file = "summary_Gini_300.txt", sep = ";", quote = FALSE, row.names = F)
```

5. P-Values 300-Year-Intervals, Gini Index

```{r}
p_values_Intervals_Gini_300 = c(1:11)

for(i in 1:11){
  
  p_values_Intervals_Gini_300[i] = Perm_Gini(X = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i])],
                                             Y = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i+1])], K = 20000, B = 999)
  
}

p_values_Intervals_Gini_300
```

6. Gini Index - Plot the results.

```{r}
graph_o_300 = table_time_series_300 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = Gini_time_series_300), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_300, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = Gini_time_series_300), color = "coral2") +
  geom_text(data = table_time_series_300, mapping = aes(x = -average_date, y = Gini_time_series_300, 
                                                        label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("G"^"*")) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3600,0), breaks = c(0,-300,-600,-900,-1200,-1500,-1800,-2100,-2400,-2700,-3000,-3300,-3600),
                     labels = c(0,300,600,900,1200,1500,1800,2100,2400,2700,3000,3300,3600),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Normalized Gini Index") +
  annotate(geom = "text", x = -3330, y = 0.6, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 0.555, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2730, y = 0.6, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.58, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2130, y = 0.71, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1830, y = 0.72, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 0.67, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 0.77, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -930, y = 0.88, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 0.93, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -330, y = 0.88, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0.4, 1), xlim = c(-3600,0), clip = "off")

graph_o_300
```

7. GEM Index (lower part) - Set up table necessary for the plot and compute the GEM index (lower part) for each interval.

```{r}
GEM_0_time_series_300 = rep(x = 0, times = length(x = absolute_periods_300))

table_time_series_300_GEM_0 = data.frame(absolute_periods_300, GEM_0_time_series_300)

table_time_series_300_GEM_0["average_date"] = c(3450,3150,2850,2550,2250,1950,1650,1350,1050,750,450,150)

table_time_series_300_GEM_0["CI_95_GEM_0_l"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300_GEM_0["CI_95_GEM_0_u"] = rep(1, times = length(x = absolute_periods_300)) 
table_time_series_300_GEM_0["GEM_0_BSE"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300_GEM_0$n = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_0[i,7] = c(length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])]))
  
}

table_time_series_300_GEM_0["sorting"] = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_0[i,2] = GEM(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])], c = 0)
  
}
```

8. GEM Index (lower part) - Compute the bootstrapped standard errors of the GEM index (lower part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_0_3600_3300 = rep(x = 0, times = B)
GEM_0_3300_3000 = rep(x = 0, times = B)
GEM_0_3000_2700 = rep(x = 0, times = B)
GEM_0_2700_2400 = rep(x = 0, times = B)
GEM_0_2400_2100 = rep(x = 0, times = B)
GEM_0_2100_1800 = rep(x = 0, times = B)
GEM_0_1800_1500 = rep(x = 0, times = B)
GEM_0_1500_1200 = rep(x = 0, times = B)
GEM_0_1200_900 = rep(x = 0, times = B)
GEM_0_900_600 = rep(x = 0, times = B)
GEM_0_600_300 = rep(x = 0, times = B)
GEM_0_300_0 = rep(x = 0, times = B)

s_1 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]
s_2 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]
s_3 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]
s_4 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]
s_5 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]
s_6 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]
s_7 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]
s_8 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]
s_9 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]
s_10 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]
s_11 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]
s_12 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3600_3300 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3300_3000 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3000_2700 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_2700_2400 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2400_2100 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2100_1800 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_1800_1500 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_1500_1200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1200_900 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_900_600 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_600_300 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_300_0 = sample(x = s_12, size = length(x = s_12), replace = T)
  
  GEM_0_3600_3300[i] = GEM(x = Sample_3600_3300, c = 0)
  GEM_0_3300_3000[i] = GEM(x = Sample_3300_3000, c = 0)
  GEM_0_3000_2700[i] = GEM(x = Sample_3000_2700, c = 0)
  GEM_0_2700_2400[i] = GEM(x = Sample_2700_2400, c = 0)
  GEM_0_2400_2100[i] = GEM(x = Sample_2400_2100, c = 0)
  GEM_0_2100_1800[i] = GEM(x = Sample_2100_1800, c = 0)
  GEM_0_1800_1500[i] = GEM(x = Sample_1800_1500, c = 0)
  GEM_0_1500_1200[i] = GEM(x = Sample_1500_1200, c = 0)
  GEM_0_1200_900[i] = GEM(x = Sample_1200_900, c = 0)
  GEM_0_900_600[i] = GEM(x = Sample_900_600, c = 0)
  GEM_0_600_300[i] = GEM(x = Sample_600_300, c = 0)
  GEM_0_300_0[i] = GEM(x = Sample_300_0, c = 0)
  
}

GEM_0_Boot_time_all_300 = rbind(GEM_0_3600_3300,GEM_0_3300_3000,GEM_0_3000_2700,GEM_0_2700_2400,GEM_0_2400_2100,GEM_0_2100_1800,GEM_0_1800_1500,
                                GEM_0_1500_1200,GEM_0_1200_900,GEM_0_900_600,GEM_0_600_300,GEM_0_300_0)

for (i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_0[i,6] = sd(x = GEM_0_Boot_time_all_300[i,])
  
}
```

9. GEM Index (lower part) - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_0[i,4] = table_time_series_300_GEM_0[i,2]-
    abs(x = quantile(x = GEM_0_Boot_time_all_300[i,], probs = c(0.025))[1])*table_time_series_300_GEM_0[i,6]
  table_time_series_300_GEM_0[i,5] = table_time_series_300_GEM_0[i,2]+
    abs(x = quantile(x = GEM_0_Boot_time_all_300[i,], probs = c(0.975))[1])*table_time_series_300_GEM_0[i,6]
  
}

table_time_series_300_GEM_0

xtable(table_time_series_300_GEM_0, type = "latex") # for Latex code

write.table(table_time_series_300_GEM_0, file = "summary_I_0_300.txt", sep = ";", quote = FALSE, row.names = F)
```

10. P-Values 300-Year-Intervals, GEM Index (lower part)

```{r}
p_values_Intervals_GEM_0_300 = c(1:11)

for(i in 1:11){
  
  p_values_Intervals_GEM_0_300[i] = Perm_GEM(X = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i])],
                                             Y = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i+1])], K = 20000, B = 999, c = 0)
  
}

p_values_Intervals_GEM_0_300
```

11. GEM Index (lower part) - Plot the results.

```{r}
graph_p_300 = table_time_series_300_GEM_0 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_0_time_series_300), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_300_GEM_0, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = GEM_0_time_series_300), color = "coral2") +
  geom_text(data = table_time_series_300_GEM_0, mapping = aes(x = -average_date, y = GEM_0_time_series_300, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[0])) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3600,0), breaks = c(0,-300,-600,-900,-1200,-1500,-1800,-2100,-2400,-2700,-3000,-3300,-3600),
                     labels = c(0,300,600,900,1200,1500,1800,2100,2400,2700,3000,3300,3600),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Generalized Entropy Measure (lower part)") +
  annotate(geom = "text", x = -3330, y = 1.15, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 0.72, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2730, y = 0.75, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.71, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2130, y = 1.14, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1830, y = 1.22, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 1.0, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 1.45, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -930, y = 2.05, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 2.36, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -330, y = 2.0, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0, 2.7), xlim = c(-3600,0), clip = "off")

graph_p_300
```

12. GEM Index (central part) - Set up table necessary for the plot and compute the GEM index (central part) for each interval.

```{r}
GEM_1_time_series_300 = rep(x = 0, times = length(x = absolute_periods_300))

table_time_series_300_GEM_1 = data.frame(absolute_periods_300, GEM_1_time_series_300)

table_time_series_300_GEM_1["average_date"] = c(3450,3150,2850,2550,2250,1950,1650,1350,1050,750,450,150)

table_time_series_300_GEM_1["CI_95_GEM_1_l"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300_GEM_1["CI_95_GEM_1_u"] = rep(1, times = length(x = absolute_periods_300)) 
table_time_series_300_GEM_1["GEM_1_BSE"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300_GEM_1$n = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_1[i,7] = c(length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])]))
  
}

table_time_series_300_GEM_1["sorting"] = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_1[i,2] = GEM(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])], c = 1)
  
}
```

13. GEM Index (central part) - Compute the bootstrapped standard errors of the GEM index (central part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_1_3600_3300 = rep(x = 0, times = B)
GEM_1_3300_3000 = rep(x = 0, times = B)
GEM_1_3000_2700 = rep(x = 0, times = B)
GEM_1_2700_2400 = rep(x = 0, times = B)
GEM_1_2400_2100 = rep(x = 0, times = B)
GEM_1_2100_1800 = rep(x = 0, times = B)
GEM_1_1800_1500 = rep(x = 0, times = B)
GEM_1_1500_1200 = rep(x = 0, times = B)
GEM_1_1200_900 = rep(x = 0, times = B)
GEM_1_900_600 = rep(x = 0, times = B)
GEM_1_600_300 = rep(x = 0, times = B)
GEM_1_300_0 = rep(x = 0, times = B)

s_1 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]
s_2 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]
s_3 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]
s_4 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]
s_5 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]
s_6 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]
s_7 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]
s_8 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]
s_9 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]
s_10 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]
s_11 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]
s_12 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3600_3300 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3300_3000 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3000_2700 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_2700_2400 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2400_2100 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2100_1800 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_1800_1500 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_1500_1200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1200_900 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_900_600 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_600_300 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_300_0 = sample(x = s_12, size = length(x = s_12), replace = T)
  
  GEM_1_3600_3300[i] = GEM(x = Sample_3600_3300, c = 1)
  GEM_1_3300_3000[i] = GEM(x = Sample_3300_3000, c = 1)
  GEM_1_3000_2700[i] = GEM(x = Sample_3000_2700, c = 1)
  GEM_1_2700_2400[i] = GEM(x = Sample_2700_2400, c = 1)
  GEM_1_2400_2100[i] = GEM(x = Sample_2400_2100, c = 1)
  GEM_1_2100_1800[i] = GEM(x = Sample_2100_1800, c = 1)
  GEM_1_1800_1500[i] = GEM(x = Sample_1800_1500, c = 1)
  GEM_1_1500_1200[i] = GEM(x = Sample_1500_1200, c = 1)
  GEM_1_1200_900[i] = GEM(x = Sample_1200_900, c = 1)
  GEM_1_900_600[i] = GEM(x = Sample_900_600, c = 1)
  GEM_1_600_300[i] = GEM(x = Sample_600_300, c = 1)
  GEM_1_300_0[i] = GEM(x = Sample_300_0, c = 1)
  
}

GEM_1_Boot_time_all_300 = rbind(GEM_1_3600_3300,GEM_1_3300_3000,GEM_1_3000_2700,GEM_1_2700_2400,GEM_1_2400_2100,GEM_1_2100_1800,GEM_1_1800_1500,
                                GEM_1_1500_1200,GEM_1_1200_900,GEM_1_900_600,GEM_1_600_300,GEM_1_300_0)

for (i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_1[i,6] = sd(x = GEM_1_Boot_time_all_300[i,])
  
}
```

14. GEM Index (central part) - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_1[i,4] = table_time_series_300_GEM_1[i,2]-
    abs(x = quantile(x = GEM_1_Boot_time_all_300[i,], probs = c(0.025))[1])*table_time_series_300_GEM_1[i,6]
  table_time_series_300_GEM_1[i,5] = table_time_series_300_GEM_1[i,2]+
    abs(x = quantile(x = GEM_1_Boot_time_all_300[i,], probs = c(0.975))[1])*table_time_series_300_GEM_1[i,6]
  
}

table_time_series_300_GEM_1

xtable(table_time_series_300_GEM_1, type = "latex") # for Latex code

write.table(table_time_series_300_GEM_1, file = "summary_I_1_300.txt", sep = ";", quote = FALSE, row.names = F)
```

15. P-Values 300-Year-Intervals, GEM Index (central part)

```{r}
p_values_Intervals_GEM_1_300 = c(1:11)

for(i in 1:11){
  
  p_values_Intervals_GEM_1_300[i] = Perm_GEM(X = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i])],
                                             Y = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i+1])], K = 20000, B = 999, c = 1)
  
}

p_values_Intervals_GEM_1_300
```

16. GEM Index (central part) - Plot the confidence intervals.

```{r}
graph_q_300 = table_time_series_300_GEM_1 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_1_time_series_300), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_300_GEM_1, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = GEM_1_time_series_300), color = "coral2") +
  geom_text(data = table_time_series_300_GEM_1, mapping = aes(x = -average_date, y = GEM_1_time_series_300, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[1])) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3600,0), breaks = c(0,-300,-600,-900,-1200,-1500,-1800,-2100,-2400,-2700,-3000,-3300,-3600),
                     labels = c(0,300,600,900,1200,1500,1800,2100,2400,2700,3000,3300,3600),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Generalized Entropy Measure (central part)") +
  annotate(geom = "text", x = -3330, y = 0.74, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 0.8, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2730, y = 0.91, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.82, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2130, y = 1.15, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1830, y = 1.15, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 1.0, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 1.65, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -930, y = 2.5, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 3.43, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -330, y = 3, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0, 3.9), xlim = c(-3600,0), clip = "off")

graph_q_300
```

17. GEM Index (upper part) - Set up table necessary for the plot and compute the GEM index (upper part) for each interval.

```{r}
GEM_2_time_series_300 = rep(x = 0, times = length(x = absolute_periods_300))

table_time_series_300_GEM_2 = data.frame(absolute_periods_300, GEM_2_time_series_300)

table_time_series_300_GEM_2["average_date"] = c(3450,3150,2850,2550,2250,1950,1650,1350,1050,750,450,150)

table_time_series_300_GEM_2["CI_95_GEM_2_l"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300_GEM_2["CI_95_GEM_2_u"] = rep(1, times = length(x = absolute_periods_300)) 
table_time_series_300_GEM_2["GEM_2_BSE"] = rep(1, times = length(x = absolute_periods_300))
table_time_series_300_GEM_2$n = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_2[i,7] = c(length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])]))
  
}

table_time_series_300_GEM_2["sorting"] = c(1:length(x = absolute_periods_300))

for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_2[i,2] = GEM(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[i])], c = 2)
  
}
```

18. GEM Index (upper part) - Compute the bootstrapped standard errors of the GEM index (upper part).

```{r}
B = 1000 # number of bootstrap samples

# Construct empty vectors for the bootstrapped inequality indices

GEM_2_3600_3300 = rep(x = 0, times = B)
GEM_2_3300_3000 = rep(x = 0, times = B)
GEM_2_3000_2700 = rep(x = 0, times = B)
GEM_2_2700_2400 = rep(x = 0, times = B)
GEM_2_2400_2100 = rep(x = 0, times = B)
GEM_2_2100_1800 = rep(x = 0, times = B)
GEM_2_1800_1500 = rep(x = 0, times = B)
GEM_2_1500_1200 = rep(x = 0, times = B)
GEM_2_1200_900 = rep(x = 0, times = B)
GEM_2_900_600 = rep(x = 0, times = B)
GEM_2_600_300 = rep(x = 0, times = B)
GEM_2_300_0 = rep(x = 0, times = B)

s_1 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]
s_2 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]
s_3 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]
s_4 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]
s_5 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]
s_6 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]
s_7 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]
s_8 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]
s_9 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]
s_10 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]
s_11 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]
s_12 = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]

set.seed(1) # since there is randomness in the draws

for(i in 1:B){
  
  Sample_3600_3300 = sample(x = s_1, size = length(x = s_1), replace = T)
  Sample_3300_3000 = sample(x = s_2, size = length(x = s_2), replace = T)
  Sample_3000_2700 = sample(x = s_3, size = length(x = s_3), replace = T)
  Sample_2700_2400 = sample(x = s_4, size = length(x = s_4), replace = T)
  Sample_2400_2100 = sample(x = s_5, size = length(x = s_5), replace = T)
  Sample_2100_1800 = sample(x = s_6, size = length(x = s_6), replace = T)
  Sample_1800_1500 = sample(x = s_7, size = length(x = s_7), replace = T)
  Sample_1500_1200 = sample(x = s_8, size = length(x = s_8), replace = T)
  Sample_1200_900 = sample(x = s_9, size = length(x = s_9), replace = T)
  Sample_900_600 = sample(x = s_10, size = length(x = s_10), replace = T)
  Sample_600_300 = sample(x = s_11, size = length(x = s_11), replace = T)
  Sample_300_0 = sample(x = s_12, size = length(x = s_12), replace = T)
  
  GEM_2_3600_3300[i] = GEM(x = Sample_3600_3300, c = 2)
  GEM_2_3300_3000[i] = GEM(x = Sample_3300_3000, c = 2)
  GEM_2_3000_2700[i] = GEM(x = Sample_3000_2700, c = 2)
  GEM_2_2700_2400[i] = GEM(x = Sample_2700_2400, c = 2)
  GEM_2_2400_2100[i] = GEM(x = Sample_2400_2100, c = 2)
  GEM_2_2100_1800[i] = GEM(x = Sample_2100_1800, c = 2)
  GEM_2_1800_1500[i] = GEM(x = Sample_1800_1500, c = 2)
  GEM_2_1500_1200[i] = GEM(x = Sample_1500_1200, c = 2)
  GEM_2_1200_900[i] = GEM(x = Sample_1200_900, c = 2)
  GEM_2_900_600[i] = GEM(x = Sample_900_600, c = 2)
  GEM_2_600_300[i] = GEM(x = Sample_600_300, c = 2)
  GEM_2_300_0[i] = GEM(x = Sample_300_0, c = 2)
  
}

GEM_2_Boot_time_all_300 = rbind(GEM_2_3600_3300,GEM_2_3300_3000,GEM_2_3000_2700,GEM_2_2700_2400,GEM_2_2400_2100,GEM_2_2100_1800,GEM_2_1800_1500,
                                GEM_2_1500_1200,GEM_2_1200_900,GEM_2_900_600,GEM_2_600_300,GEM_2_300_0)

for (i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_2[i,6] = sd(x = GEM_2_Boot_time_all_300[i,])
  
}
```

19. GEM Index (upper part) - Compute the confidence intervals.

```{r}
for(i in 1:length(x = absolute_periods_300)){
  
  table_time_series_300_GEM_2[i,4] = table_time_series_300_GEM_2[i,2]-
    abs(x = quantile(x = GEM_2_Boot_time_all_300[i,], probs = c(0.025))[1])*table_time_series_300_GEM_2[i,6]
  table_time_series_300_GEM_2[i,5] = table_time_series_300_GEM_2[i,2]+
    abs(x = quantile(x = GEM_2_Boot_time_all_300[i,], probs = c(0.975))[1])*table_time_series_300_GEM_2[i,6]
  
}

table_time_series_300_GEM_2

xtable(table_time_series_300_GEM_2, type = "latex") # for Latex code

write.table(table_time_series_300_GEM_2, file = "summary_I_2_300.txt", sep = ";", quote = FALSE, row.names = F)
```

20. P-Values 300-Year-Intervals, GEM Index (upper part)

```{r}
p_values_Intervals_GEM_2_300 = c(1:11)

for(i in 1:11){
  
  p_values_Intervals_GEM_2_300[i] = Perm_GEM(X = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i])],
                                             Y = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in%                                                 absolute_periods_300[i+1])], K = 20000, B = 999, c = 2)
  
}

p_values_Intervals_GEM_2_300
```

21. GEM Index (upper part) - Plot the results.

```{r}
graph_r_300 = table_time_series_300_GEM_2 %>%
  
  ggplot() + 
  geom_point(mapping = aes(x = -average_date, y = GEM_2_time_series_300), size = 4.5, color = "coral2") + 
  geom_line(data = table_time_series_300_GEM_2, size = 2, linetype = "solid",
            mapping = aes(x = -average_date, y = GEM_2_time_series_300), color = "coral2") +
  geom_text(data = table_time_series_300_GEM_2, mapping = aes(x = -average_date, y = GEM_2_time_series_300, 
                                                              label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 10, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 15, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank()) +
  ylab(expression("I"[2])) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3600,0), breaks = c(0,-300,-600,-900,-1200,-1500,-1800,-2100,-2400,-2700,-3000,-3300,-3600),
                     labels = c(0,300,600,900,1200,1500,1800,2100,2400,2700,3000,3300,3600),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  ggtitle("Generalized Entropy Measure (upper part)") +
  annotate(geom = "text", x = -3330, y = 3.5, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -3030, y = 4.8, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2730, y = 5.3, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 4.7, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2130, y = 5.3, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1830, y = 5.1, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1530, y = 4.8, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1230, y = 10.5, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -930, y = 23, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 61, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -330, y = 54, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  coord_cartesian(ylim = c(0, 80), xlim = c(-3600,0), clip = "off")

graph_r_300
```

22. Plot all graphs in one figure.

```{r}
pdf(file = "path.pdf",
    width = 34, height = 16)

ggarrange(graph_o_300, graph_p_300, graph_q_300, graph_r_300, ncol = 2, nrow = 2)

dev.off()
```

23. Produce descriptive statistics for every 300-year interval.

```{r}
# Set up data frame for summary statistics

stats_3600_3300 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[1])]))

stats_3300_3000 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[2])]))

stats_3000_2700 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[3])]))

stats_2700_2400 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[4])]))

stats_2400_2100 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[5])]))

stats_2100_1800 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[6])]))

stats_1800_1500 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[7])]))

stats_1500_1200 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[8])]))

stats_1200_900 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]),
                    max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]),
                    quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])],
                             probs = c(0.25,0.5,0.75)),
                    mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]),
                    sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]),
                    length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[9])]))

stats_900_600 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]),
                   max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]),
                   quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])],
                            probs = c(0.25,0.5,0.75)),
                   mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]),
                   sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]),
                   length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[10])]))

stats_600_300 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]),
                  max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]),
                  quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])],
                           probs = c(0.25,0.5,0.75)),
                  mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]),
                  sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]),
                  length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[11])]))

stats_300_0 = c(min(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]),
                  max(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]),
                  quantile(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])],
                           probs = c(0.25,0.5,0.75)),
                  mean(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]),
                  sd(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]),
                  length(x = data_time_300$barrow_volume[which(data_time_300[,"time_series_dating"] %in% absolute_periods_300[12])]))


summary_statistics_300_years = round(x =data.frame(rbind(stats_3600_3300,stats_3300_3000,stats_3000_2700,stats_2700_2400,stats_2400_2100,
                                                         stats_2100_1800,stats_1800_1500,stats_1500_1200,stats_1200_900,stats_900_600,
                                                         stats_600_300,stats_300_0)),digit = 2)

colnames(summary_statistics_300_years) = c("min","max","q_0.25","q_0.5","q_0.75","mean","std_dev","num_obs")

xtable(summary_statistics_300_years, type = "latex") # for Latex code

write.table(summary_statistics_300_years, file = "summary_300.txt", sep = ";", quote = FALSE, row.names = F)
```



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



Section 4: Main Figure of the Analysis



####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 
####################################################################################################### 



1. Import the dataset for the share of burial mounds chart.

```{r}
dataset_share_burial_mounds = read_excel("dataset_share_burial_mounds.xlsx")

dataset_share_burial_mounds_larger_1400 = read_excel("dataset_share_burial_mounds_larger_1400.xlsx") # for the appendix graph
```

2. Joint plot of the Gini Index, total burials, and share of burial mounds in total burials (main figure).

```{r}
scale_factor_sum_burials = 1/max(x = dataset_share_burial_mounds$number_buried_sum)

pdf(file = "path.pdf",
    width = 34, height = 16)

joint_plot_ineq =  ggplot() + 
  geom_bar(data = dataset_share_burial_mounds, mapping = aes(x = BCE, y = number_buried_sum*scale_factor_sum_burials,
                                                             fill = "Sum of Buried Individuals"),
           stat = "identity", color = NA, alpha = 0.65) +
  geom_point(data = table_time_series_200, mapping = aes(x = -average_date, y = Gini_time_series_200), size = 4.5,
             color = "coral2") + 
  geom_point(data = dataset_share_burial_mounds,
             mapping = aes(x = BCE, y = share_burial_mound_buried), size = 4.5, color = "black") + 
  geom_line(data = table_time_series_200, size = 2,
            mapping = aes(x = -average_date, y = Gini_time_series_200, colour = "Gini Index"), linetype = "solid") +
  geom_line(data = dataset_share_burial_mounds, size = 2,
            mapping = aes(x = BCE, y = share_burial_mound_buried,
                          colour = "Share Buried in Burial Mounds (Relative Size of the Upper Societal Segment)"),
            linetype = "solid") +
  geom_text(data = table_time_series_200, mapping = aes(x = -average_date, y = Gini_time_series_200, 
                                                        label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 50, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 50, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        axis.title.y.right = element_text(colour = "black", margin = margin(t = 0, r = 70, b = 0, l = 0), vjust = 10,
                                          size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank(),
        legend.title = element_text(size = 36),
        legend.text = element_text(size = 30),
        legend.spacing.y = unit(0.5, "cm"),
        legend.key.size = unit(3,"line"),
        legend.key.height = unit(4.5,"line"),
        legend.position = "bottom") +
  ylab(expression("Gini Index and Share Buried in Burial Mounds")) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3800,-200),
                     breaks = c(-200,-400,-600,-800,-1000,-1200,-1400,-1600,-1800,-2000,-2200,-2400,-2600,-2800,-3000,-3200,
                                -3400,-3600,-3800),
                     labels = c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(limits = c(0,1.08),
                     breaks = seq(from = 0, to = 1, by = 0.1),
                     labels = seq(from = 0, to = 1, by = 0.1),
                     sec.axis = sec_axis(~./scale_factor_sum_burials,
                                         labels = c("0","5,000","10,000","15,000","20,000","25,000","30,000","35,000",
                                                    "40,000","45,000"),
                                         breaks = seq(from = 0, to = 45000, by = 5000),
                                         name = "Sum of Buried Individuals")) +
  scale_colour_manual("", breaks = c("Gini Index",
                                     "Share Buried in Burial Mounds (Relative Size of the Upper Societal Segment)",
                                     "Sum of Buried Individuals"),
                      values = c("Gini Index"="coral2","Share Buried in Burial Mounds (Relative Size of the Upper Societal Segment)"="black")) +
  #scale_fill_manual(name = "", values = c("Sum of Buried Individuals"="forestgreen")) +
  annotate(geom = "text", x = -450, y = 0.905, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -630, y = 0.885, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -850, y = 0.81, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1030, y = 0.795, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1250, y = 0.775, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1430, y = 0.69, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -1650, y = 0.64, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1830, y = 0.79, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2030, y = 0.77, label = "x", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2250, y = 0.575, label = "o", hjust = "left", size = 12, colour = "coral2") + 
  annotate(geom = "text", x = -2430, y = 0.51, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2630, y = 0.588, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3030, y = 0.50, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3430, y = 0.64, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3630, y = 0.70, label = "o", hjust = "left", size = 12, colour = "coral2") +
  geom_rect(
    aes(xmin = -3800, xmax = -3500, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -3500, xmax = -3300, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -3300, xmax = -2800, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -2800, xmax = -2200, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -2200, xmax = -1600, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -1600, xmax = -1200, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -1200, xmax = -800, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -800, xmax = -400, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -400, xmax = -200, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  scale_fill_manual(
    values = c("Phase_1" = "floralwhite", "Phase_2" = "grey93","Sum of Buried Individuals"="forestgreen"),
    guide = "none"
  ) +
  annotate("rect", xmin = -3800, xmax = -3500, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -3650, y = 1.05, label = "Ph. 1.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -3500, xmax = -3300, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -3400, y = 1.05, label = "Ph. 1.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -3300, xmax = -2800, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -3050, y = 1.05, label = "Ph. 2.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -2800, xmax = -2200, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -2500, y = 1.05, label = "Ph. 2.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -2200, xmax = -1600, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -1900, y = 1.05, label = "Ph. 3.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -1600, xmax = -1200, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -1400, y = 1.05, label = "Ph. 3.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -1200, xmax = -800, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -1000, y = 1.05, label = "Ph. 4.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -800, xmax = -400, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -600, y = 1.05, label = "Ph. 4.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -400, xmax = -200, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -300, y = 1.05, label = "...", hjust = 0.5, size = 10) +
  geom_segment(aes(x = -3800, xend = -3800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -3800, xend = -3800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -3500, xend = -3500, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -3300, xend = -3300, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -2800, xend = -2800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -2200, xend = -2200, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -1600, xend = -1600, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -1200, xend = -1200, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -800, xend = -800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -400, xend = -400, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1)

joint_plot_ineq

dev.off()
```

3. Joint plot of the Gini Index, total burials, and share of burial mounds in total burials (1,400 years, Appendix; only works properly if you set the exclusion restriction in Section 2.1 to >1400).

```{r}
scale_factor_sum_burials = 1/max(x = dataset_share_burial_mounds_larger_1400$number_buried_sum)

pdf(file = "path.pdf",
    width = 34, height = 16)

joint_plot_ineq_1400 =  ggplot() +
  geom_bar(data = dataset_share_burial_mounds_larger_1400, mapping = aes(x = BCE, y = number_buried_sum*scale_factor_sum_burials,
                                                             fill = "Sum of Buried Individuals"),
           stat = "identity", color = NA, alpha = 0.65) +
  geom_point(data = table_time_series_200, mapping = aes(x = -average_date, y = Gini_time_series_200), size = 4.5,
             color = "coral2") +
  geom_point(data = dataset_share_burial_mounds_larger_1400,
             mapping = aes(x = BCE, y = share_burial_mound_buried), size = 4.5, color = "black") +
  geom_line(data = table_time_series_200, size = 2,
            mapping = aes(x = -average_date, y = Gini_time_series_200, colour = "Gini Index"), linetype = "solid") +
  geom_line(data = dataset_share_burial_mounds_larger_1400, size = 2,
            mapping = aes(x = BCE, y = share_burial_mound_buried,
                          colour = "Share Buried in Burial Mounds (Relative Size of the Upper Societal Segment)"),
            linetype = "solid") +
  geom_text(data = table_time_series_200, mapping = aes(x = -average_date, y = Gini_time_series_200,
                                                        label = n), vjust = -0.6, size = 9) +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5, size = 36),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1, size = 28),
        axis.title.x = element_text(colour = "black", margin = margin(t = 50, r = 0, b = 0, l = 0), size = 28),
        axis.title.y = element_text(colour = "black", margin = margin(t = 0, r = 50, b = 0, l = 0), size = 28),
        axis.text.y = element_text(angle = 0, vjust = 0.5, hjust = 1, size = 28),
        axis.text = element_text(color = "black", size = 28),
        axis.title.y.right = element_text(colour = "black", margin = margin(t = 0, r = 70, b = 0, l = 0), vjust = 10,
                                          size = 28),
        panel.grid.major = element_line(colour = "grey95", size = 0.2),
        panel.grid.minor = element_blank(),
        legend.title = element_text(size = 36),
        legend.text = element_text(size = 30),
        legend.spacing.y = unit(0.5, "cm"),
        legend.key.size = unit(3,"line"),
        legend.key.height = unit(4.5,"line"),
        legend.position = "bottom") +
  ylab(expression("Gini Index and Share Buried in Burial Mounds")) +
  xlab("Years BCE") +
  scale_x_continuous(limits = c(-3800,-200),
                     breaks = c(-200,-400,-600,-800,-1000,-1200,-1400,-1600,-1800,-2000,-2200,-2400,-2600,-2800,-3000,-3200,
                                -3400,-3600,-3800),
                     labels = c(200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,2800,3000,3200,3400,3600,3800),
                     sec.axis = sec_axis(~., labels = NULL, breaks = NULL)) +
  scale_y_continuous(limits = c(0,1.08),
                     breaks = seq(from = 0, to = 1, by = 0.1),
                     labels = seq(from = 0, to = 1, by = 0.1),
                     sec.axis = sec_axis(~./scale_factor_sum_burials,
                                         labels = c("0","5,000","10,000","15,000","20,000","25,000","30,000","35,000",
                                                    "40,000","45,000"),
                                         breaks = seq(from = 0, to = 45000, by = 5000),
                                         name = "Sum of Buried Individuals")) +
  scale_colour_manual("", breaks = c("Gini Index",
                                     "Share Buried in Burial Mounds (Relative Size of the Upper Societal Segment)",
                                     "Sum of Buried Individuals"),
                      values = c("Gini Index"="coral2","Share Buried in Burial Mounds (Relative Size of the Upper Societal Segment)"="black")) +
  #scale_fill_manual(name = "", values = c("Sum of Buried Individuals"="forestgreen")) +
  annotate(geom = "text", x = -430, y = 0.93, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -630, y = 0.89, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -830, y = 0.89, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1030, y = 0.8, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1230, y = 0.735, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1430, y = 0.69, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1630, y = 0.65, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -1830, y = 0.79, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2030, y = 0.78, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2230, y = 0.58, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2430, y = 0.51, label = "o", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -2630, y = 0.588, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3030, y = 0.50, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3430, y = 0.57, label = "x", hjust = "left", size = 12, colour = "coral2") +
  annotate(geom = "text", x = -3630, y = 0.67, label = "o", hjust = "left", size = 12, colour = "coral2") +
  geom_rect(
    aes(xmin = -3800, xmax = -3500, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -3500, xmax = -3300, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -3300, xmax = -2800, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -2800, xmax = -2200, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -2200, xmax = -1600, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -1600, xmax = -1200, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -1200, xmax = -800, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  geom_rect(
    aes(xmin = -800, xmax = -400, ymin = 0, ymax = 1.03, fill = "Phase_2"),
    alpha = 0.25
  ) +
  geom_rect(
    aes(xmin = -400, xmax = -200, ymin = 0, ymax = 1.03, fill = "Phase_1"),
    alpha = 0.2
  ) +
  scale_fill_manual(
    values = c("Phase_1" = "floralwhite", "Phase_2" = "grey93","Sum of Buried Individuals"="forestgreen"),
    guide = "none"
  ) +
  annotate("rect", xmin = -3800, xmax = -3500, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -3650, y = 1.05, label = "Ph. 1.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -3500, xmax = -3300, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -3400, y = 1.05, label = "Ph. 1.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -3300, xmax = -2800, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -3050, y = 1.05, label = "Ph. 2.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -2800, xmax = -2200, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -2500, y = 1.05, label = "Ph. 2.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -2200, xmax = -1600, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -1900, y = 1.05, label = "Ph. 3.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -1600, xmax = -1200, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -1400, y = 1.05, label = "Ph. 3.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -1200, xmax = -800, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -1000, y = 1.05, label = "Ph. 4.1", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -800, xmax = -400, ymin = 1.03, ymax = 1.07, fill = "grey93", alpha = 1) +
  annotate("text", x = -600, y = 1.05, label = "Ph. 4.2", hjust = 0.5, size = 10) +
  annotate("rect", xmin = -400, xmax = -200, ymin = 1.03, ymax = 1.07, fill = "floralwhite", alpha = 1) +
  annotate("text", x = -300, y = 1.05, label = "...", hjust = 0.5, size = 10) +
  geom_segment(aes(x = -3800, xend = -3800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -3800, xend = -3800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -3500, xend = -3500, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -3300, xend = -3300, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -2800, xend = -2800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -2200, xend = -2200, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -1600, xend = -1600, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -1200, xend = -1200, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -800, xend = -800, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1) +
  geom_segment(aes(x = -400, xend = -400, y = 0, yend = 1.07), linetype = "dashed", alpha = 0.5,
               col = "black", size = 1.1)

joint_plot_ineq_1400

dev.off()
```
